!************************
MODULE libinterp
!***********************
! User definitions of structures
USE ElliDef, ONLY : Precision

! USE problem ! the external routines used to calculate coefficients
             ! in the poisson equations are stored here

 IMPLICIT NONE
 SAVE
 PRIVATE
 PUBLIC rgrd1,rgrd2,rgrd3
INTEGER, PARAMETER :: qPrec=Precision
! 
 CONTAINS

   !
 
! Code converted using TO_F90 by Alan Miller
! Date: 2006-11-29  Time: 13:21:04

! ... file rgrd1.f

!     this file contains documentation for subroutine rgrd1 followed by
!     fortran code for rgrd1 and additional subroutines.

! ... author

!     John C. Adams (NCAR 1997)

! ... subroutine rgrd1(nx,x,p,mx,xx,q,intpol,w,lw,iw,liw,ier)

! ... purpose

!     subroutine rgrd1 interpolates the values p(i) on the grid x(i)
!     for i=1,...,nx onto q(ii) on the grid xx(ii),ii=1,...,mx.

! ... language

!     coded in portable FORTRAN90 and FORTRAN77

! ... test program

!     file trgrd1.f on regridpack includes a test program for subroutine rgrd1

! ... method

!     linear or cubic interpolation is used (see argument intpol)

! ... required files

!     none

! ... requirements

!     x must be a strictly increasing grid and xx must be an increasing
!     grid (see ier = 4).  in addition the interval

!          [xx(1),xx(mx)]

!     must lie within the interval

!          [x(1),x(nx)].

!     extrapolation is not allowed (see ier=3).  if these intervals
!     are identical and the x and xx grids are UNIFORM then subroutine
!     rgrd1u (see file rgrd1u.f) should be used in place of rgrd1.

! ... required files

!     none

! ... efficiency

!     inner most loops in regridpack software vectorize.


! *** input arguments


! ... nx

!     the integer dimension of the grid vector x and the dimension of p.
!     nx > 1 if intpol = 1 or nx > 3 if intpol = 3 is required.

! ... x

!     a real nx vector of strictly increasing values which defines the x
!     grid on which p is given.


! ... p

!     a real nx vector of values given on the x grid

! ... mx

!     the integer dimension of the grid vector xx and the dimension of q.
!     mx > 0 is required.

! ... xx

!     a real mx vector of increasing values which defines the
!     grid on which q is defined.  xx(1) < x(1) or xx(mx) > x(nx)
!     is not allowed (see ier = 3)

! ... intpol

!     an integer which sets linear or cubic
!     interpolation as follows:

!        intpol = 1 sets linear interpolation
!        intpol = 3 sets cubic interpolation

!     values other than 1 or 3 in intpol are not allowed (ier = 6).

! ... w

!     a real work space of length at least lw which must be provided in the
!     routine calling rgrd1


! ... lw

!     the integer length of the real work space w.  let

!          lwmin =   mx            if intpol(1) = 1
!          lwmin = 4*mx            if intpol(1) = 3

!     then lw must be greater than or equal to lwmin

! ... iw

!     an integer work space of length at least liw which must be provided in the
!     routine calling rgrd1

! ... liw

!     tne length of the integer work space iw. liw must be greater than or equal to mx.


! *** output arguments


! ... q

!     a real mx vector of values on the xx grid which are
!     interpolated from p on the x grid

! ... ier

!     an integer error flag set as follows:

!     ier = 0 if no errors in input arguments are detected

!     ier = 1 if  mx < 1

!     ier = 2 if nx < 2 when intpol=1 or nx < 4 when intpol=3

!     ier = 3 if xx(1) < x(1) or x(nx) < xx(mx)

! *** to avoid this flag when end points are intended to be the
!     same but may differ slightly due to roundoff error, they
!     should be set exactly in the calling routine (e.g., if both
!     grids have the same x boundaries then xx(1)=x(1) and xx(mx)=x(nx)
!     should be set before calling rgrd1)

!     ier = 4 if the x grid is not strictly monotonically increasing
!             or if the xx grid is not montonically increasing.  more
!             precisely if:

!             x(i+1) <= x(i) for some i such that 1 <= i < nx (or)

!             xx(ii+1) < xx(ii) for some ii such that 1 <= ii < mx

!     ier = 5 if lw or liw is too small (insufficient work space)

!     ier = 6 if intpol is not equal to 1 or 3

! ************************************************************************

!     end of rgrd1 documentation, fortran code follows:

! ************************************************************************

SUBROUTINE rgrd1(nx,x,p,mx,xx,q,intpol,w,lw,iw,liw,ier)
!     dimension x(nx),p(nx),xx(mx),q(mx),w(lw)
IMPLICIT NONE
INTEGER, INTENT(IN OUT)                      :: nx
REAL (KIND=qPrec), INTENT(IN OUT)                         :: x(*)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: p(*)
INTEGER, INTENT(IN OUT)                      :: mx
REAL (KIND=qPrec), INTENT(IN OUT)                         :: xx(*)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: q(*)
INTEGER, INTENT(IN OUT)                      :: intpol
REAL (KIND=qPrec), INTENT(IN OUT)                     :: w(*)
INTEGER, INTENT(IN OUT)                  :: lw
INTEGER, INTENT(IN OUT)                  :: iw(*)
INTEGER, INTENT(IN OUT)                  :: liw
INTEGER, INTENT(IN OUT)                     :: ier



INTEGER :: i,ii,i1,i2,i3,i4

!     check arguments for errors

ier = 1

!     check xx grid resolution

IF (mx < 1) RETURN

!     check intpol

ier = 6
IF (intpol /= 1 .AND. intpol /= 3) RETURN

!     check x grid resolution

ier = 2
IF (intpol == 1 .AND. nx < 2) RETURN
IF (intpol == 3 .AND. nx < 4) RETURN

!     check xx grid contained in x grid

ier = 3
IF (xx(1) < x(1) .OR. xx(mx) > x(nx)) RETURN

!     check montonicity of grids

DO i=2,nx
  IF (x(i-1) >= x(i)) THEN
    ier = 4
    RETURN
  END IF
END DO
DO ii=2,mx
  IF (xx(ii-1) > xx(ii)) THEN
    ier = 4
    RETURN
  END IF
END DO

!     check minimum work space lengths

IF (intpol == 1) THEN
  IF (lw < mx) RETURN
ELSE
  IF (lw < 4*mx) RETURN
END IF
IF (liw < mx) RETURN

!     arguments o.k.

ier = 0

IF (intpol == 1) THEN
  
!     linear interpolation in x
  
  CALL linmx(nx,x,mx,xx,iw,w)
  CALL lint1(nx,p,mx,q,iw,w)
  RETURN
ELSE
  
!     cubic interpolation in x
  
  i1 = 1
  i2 = i1+mx
  i3 = i2+mx
  i4 = i3+mx
  CALL cubnmx(nx,x,mx,xx,iw,w(i1),w(i2),w(i3),w(i4))
  CALL cubt1(nx,p,mx,q,iw,w(i1),w(i2),w(i3),w(i4))
  RETURN
END IF
END SUBROUTINE rgrd1

SUBROUTINE lint1(nx,p,mx,q,ix,dx)
IMPLICIT NONE
INTEGER, INTENT(IN )                  :: nx

REAL (KIND=qPrec), INTENT(IN)                         :: p(nx)
INTEGER, INTENT(IN)                      :: mx
REAL (KIND=qPrec), INTENT(IN OUT)                        :: q(mx)
INTEGER, INTENT(IN OUT)                      :: ix(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                         :: dx(mx)

INTEGER :: ii,i


!     linearly interpolate p on x onto q on xx

DO ii=1,mx
  i = ix(ii)
  q(ii) = p(i)+dx(ii)*(p(i+1)-p(i))
END DO
RETURN
END SUBROUTINE lint1

SUBROUTINE cubt1(nx,p,mx,q,ix,dxm,dx,dxp,dxpp)
IMPLICIT NONE
INTEGER, INTENT(IN )                  :: nx
REAL (KIND=qPrec), INTENT(IN)                         :: p(nx)
INTEGER, INTENT(IN)                      :: mx
REAL (KIND=qPrec), INTENT(IN OUT)                        :: q(mx)
INTEGER, INTENT(IN OUT)                      :: ix(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                         :: dxm(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                         :: dx(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                         :: dxp(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                         :: dxpp(mx)

INTEGER :: i,ii


!     cubically interpolate p on x to q on xx

DO ii=1,mx
  i = ix(ii)
  q(ii) = dxm(ii)*p(i-1)+dx(ii)*p(i)+dxp(ii)*p(i+1)+dxpp(ii)*p(i+2)
END DO
RETURN
END SUBROUTINE cubt1

SUBROUTINE linmx(nx,x,mx,xx,ix,dx)

!     set x grid pointers for xx grid and interpolation scale terms

IMPLICIT NONE
INTEGER, INTENT(IN )                      :: nx
REAL (KIND=qPrec), INTENT(IN)                         :: x(*)
INTEGER, INTENT(IN)                      :: mx
REAL (KIND=qPrec), INTENT(IN)                         :: xx(*)
INTEGER, INTENT(IN OUT)                     :: ix(*)
REAL (KIND=qPrec), INTENT(IN OUT)                        :: dx(*)


INTEGER :: isrt,ii,i

isrt = 1
DO ii=1,mx
  
!     find x(i) s.t. x(i) < xx(ii) <= x(i+1)
  
  DO i=isrt,nx-1
    IF (x(i+1) >= xx(ii)) THEN
      isrt = i
      ix(ii) = i
      GO TO 3
    END IF
  END DO
  3   CONTINUE
END DO

!     set linear scale term

DO ii=1,mx
  i = ix(ii)
  dx(ii) = (xx(ii)-x(i))/(x(i+1)-x(i))
END DO
RETURN
END SUBROUTINE linmx

SUBROUTINE cubnmx(nx,x,mx,xx,ix,dxm,dx,dxp,dxpp)
IMPLICIT NONE
INTEGER, INTENT(IN )                      :: nx
REAL (KIND=qPrec), INTENT(IN)                         :: x(*)
INTEGER, INTENT(IN )                      :: mx
REAL (KIND=qPrec), INTENT(IN)                         :: xx(*)
INTEGER, INTENT(IN OUT)                     :: ix(*)
REAL (KIND=qPrec), INTENT(IN OUT)                        :: dxm(*)
REAL (KIND=qPrec), INTENT(IN OUT)                        :: dx(*)
REAL (KIND=qPrec), INTENT(IN OUT)                        :: dxp(*)
REAL (KIND=qPrec), INTENT(IN OUT)                        :: dxpp(*)


INTEGER :: i,ii,isrt

isrt = 1
DO ii=1,mx
  
!     set i in [2,nx-2] closest s.t.
!     x(i-1),x(i),x(i+1),x(i+2) can interpolate xx(ii)
  
  DO i=isrt,nx-1
    IF (x(i+1) >= xx(ii)) THEN
      ix(ii) = MIN0(nx-2,MAX0(2,i))
      isrt = ix(ii)
      GO TO 3
    END IF
  END DO
  3   CONTINUE
END DO

!     set cubic scale terms

DO ii=1,mx
  i = ix(ii)
  dxm(ii) = (xx(ii)-x(i))*(xx(ii)-x(i+1))*(xx(ii)-x(i+2))/  &
      ((x(i-1)-x(i))*(x(i-1)-x(i+1))*(x(i-1)-x(i+2)))
  dx(ii) = (xx(ii)-x(i-1))*(xx(ii)-x(i+1))*(xx(ii)-x(i+2))/  &
      ((x(i)-x(i-1))*(x(i)-x(i+1))*(x(i)-x(i+2)))
  dxp(ii) = (xx(ii)-x(i-1))*(xx(ii)-x(i))*(xx(ii)-x(i+2))/  &
      ((x(i+1)-x(i-1))*(x(i+1)-x(i))*(x(i+1)-x(i+2)))
  dxpp(ii) = (xx(ii)-x(i-1))*(xx(ii)-x(i))*(xx(ii)-x(i+1))/  &
      ((x(i+2)-x(i-1))*(x(i+2)-x(i))*(x(i+2)-x(i+1)))
END DO
RETURN
END SUBROUTINE cubnmx

!
 
! Code converted using TO_F90 by Alan Miller
! Date: 2006-11-29  Time: 13:21:09

! ... file rgrd2.f

!     this file contains documentation for subroutine rgrd2 followed by
!     fortran code for rgrd2 and additional subroutines.

! ... author

!     John C. Adams (NCAR 1997)

! ... subroutine rgrd2(nx,ny,x,y,p,mx,my,xx,yy,q,intpol,w,lw,iw,liw,ier)

! ... purpose

!     subroutine rgrd2 interpolates the values p(i,j) on the orthogonal
!     grid (x(i),y(j)) for i=1,...,nx and j=1,...,ny onto q(ii,jj) on the
!     orthogonal grid (xx(ii),yy(jj)) for ii=1,...,mx and jj=1,...,my.

! ... language

!     coded in portable FORTRAN77 and FORTRAN90

! ... test program

!     file trgrd2.f on regridpack includes a test program for subroutine rgrd2

! ... method

!     linear or cubic interpolation is used (independently) in
!     each direction (see argument intpol).

! ... required files

!     file rgrd1.f must be loaded with rgrd2.f.  it includes
!     subroutines called by the routines in rgrd2.f

! ... requirements

!     each of the x,y grids must be strictly montonically increasing
!     and each of the xx,yy grids must be montonically increasing (see
!     ier = 4).  in addition the (X,Y) region

!          [xx(1),xx(mx)] X [yy(1),yy(my)]

!     must lie within the (X,Y) region

!          [x(1),x(nx)] X [y(1),y(ny)].

!     extrapolation is not allowed (see ier=3).  if these (X,Y)
!     regions are identical and the orthogonal grids are UNIFORM
!     in each direction then subroutine rgrd2u (see file rgrd2u.f)
!     should be used instead of rgrd2.

! ... efficiency

!     inner most loops in regridpack software vectorize.
!     If the arguments mx,my (see below) have different values, optimal
!     vectorization will be achieved if mx > my.


! *** input argument


! ... nx

!     the integer dimension of the grid vector x and the first dimension
!     of p.  nx > 1 if intpol(1) = 1 or nx > 3 if intpol(1) = 3 is required.

! ... ny

!     the integer dimension of the grid vector y and the second dimension
!     of p.  ny > 1 if intpol(2) = 1 or ny > 3 if intpol(2) = 3 is required.

! ... x

!     a real nx vector of strictly increasing values which defines the x
!     portion of the orthogonal grid on which p is given

! ... y

!     a real ny vector of strictly increasing values which defines the y
!     portion of the orthogonal grid on which p is given

! ... p

!     a real nx by ny array of values given on the orthogonal (x,y) grid

! ... mx

!     the integer dimension of the grid vector xx and the first dimension
!     of q.  mx > 0 is required.

! ... my

!     the integer dimension of the grid vector yy and the second dimension
!     of q.  my > 0 is required.

! ... xx

!     a real mx vector of increasing values which defines the x portion of the
!     orthogonal grid on which q is defined.  xx(1) < x(1) or xx(mx) > x(nx)
!     is not allowed (see ier = 3)

! ... yy

!     a real my vector of increasing values which defines the y portion of the
!     orthogonal grid on which q is defined.  yy(1) < y(1) or yy(my) > y(ny)
!     is not allowed (see ier = 3)

! ... intpol

!     an integer vector of dimension 2 which sets linear or cubic
!     interpolation in the x,y directions as follows:

!        intpol(1) = 1 sets linear interpolation in the x direction
!        intpol(1) = 3 sets cubic interpolation in the x direction.

!        intpol(2) = 1 sets linear interpolation in the y direction
!        intpol(2) = 3 sets cubic interpolation in the y direction.

!     values other than 1 or 3 in intpol are not allowed (ier = 5).

! ... w

!     a real work space of length at least lw which must be provided in the
!     routine calling rgrd2

! ... lw

!     the integer length of the real work space w.  let

!          lwx =   mx            if intpol(1) = 1
!          lwx = 4*mx            if intpol(1) = 3

!          lwy = my+2*mx         if intpol(2) = 1
!          lwy = 4*(mx+my)       if intpol(2) = 3

!     then lw must be greater than or equal to lwx+lwy

! ... iw

!     an integer work space of length at least liw which must be provided in the
!     routine calling rgrd2

! ... liw

!     the integer length of the integer work space iw.  liw must be at least mx+my

! *** output arguments


! ... q

!     a real mx by my array of values on the (xx,yy) grid which are
!     interpolated from p on the (x,y) grid

! ... ier

!     an integer error flag set as follows:

!     ier = 0 if no errors in input arguments are detected

!     ier = 1 if  min0(mx,my) < 1

!     ier = 2 if nx < 2 when intpol(1)=1 or nx < 4 when intpol(1)=3 (or)
!                ny < 2 when intpol(2)=1 or ny < 4 when intpol(2)=3

!     ier = 3 if xx(1) < x(1) or x(nx) < xx(mx) (or)
!                yy(1) < y(1) or y(ny) < yy(my) (or)

! *** to avoid this flag when end points are intended to be the
!     same but may differ slightly due to roundoff error, they
!     should be set exactly in the calling routine (e.g., if both
!     grids have the same y boundaries then yy(1)=y(1) and yy(my)=y(ny)
!     should be set before calling rgrd2)

!     ier = 4 if one of the grids x,y is not strictly monotonically
!             increasing or if one of the grids xx,yy is not
!             montonically increasing.  more precisely if:

!             x(i+1) <= x(i) for some i such that 1 <= i < nx (or)

!             y(j+1) <= y(j) for some j such that 1 <= j < ny (or)

!             xx(ii+1) < xx(ii) for some ii such that 1 <= ii < mx (or)

!             yy(jj+1) < yy(jj) for some jj such that 1 <= jj < my

!     ier = 5 if lw or liw is to small (insufficient work space)

!     ier = 6 if intpol(1) or intpol(2) is not equal to 1 or 3

! ************************************************************************

!     end of rgrd2 documentation, fortran code follows:

! ************************************************************************

SUBROUTINE rgrd2(nx,ny,x,y,p,mx,my,xx,yy,q,intpol,w,lw,iw,liw,ier)
IMPLICIT NONE
INTEGER, INTENT(IN )                      :: nx
INTEGER, INTENT(IN )                      :: ny
REAL (KIND=qPrec), INTENT(IN)                         :: x(nx)
REAL (KIND=qPrec), INTENT(IN)                         :: y(ny)
REAL (KIND=qPrec), INTENT(IN)                     :: p(nx,ny)
INTEGER, INTENT(IN)                      :: mx
INTEGER, INTENT(IN)                      :: my
REAL (KIND=qPrec), INTENT(IN)                         :: xx(mx)
REAL (KIND=qPrec), INTENT(IN)                         :: yy(my)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: q(mx,my)
INTEGER, INTENT(IN)                      :: intpol(2)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: w(lw)
INTEGER, INTENT(IN)                  :: lw
INTEGER, INTENT(IN OUT)                  :: iw(liw)
INTEGER, INTENT(IN)                  :: liw
INTEGER, INTENT(IN OUT)                     :: ier




INTEGER :: i,ii,j,jj,j2,j3,j4,j5,j6,j7,j8,j9,i2,i3,i4,i5
INTEGER :: jy,lwx,lwy

!     check input arguments

ier = 1

!     check (xx,yy) grid resolution

IF (MIN0(mx,my) < 1) RETURN

!     check intpol

ier = 6
IF (intpol(1) /= 1 .AND. intpol(1) /= 3) RETURN
IF (intpol(2) /= 1 .AND. intpol(2) /= 3) RETURN

!     check (x,y) grid resolution

ier = 2
IF (intpol(1) == 1 .AND. nx < 2) RETURN
IF (intpol(1) == 3 .AND. nx < 4) RETURN
IF (intpol(2) == 1 .AND. ny < 2) RETURN
IF (intpol(2) == 3 .AND. ny < 4) RETURN

!     check work space lengths

ier = 5
IF (intpol(1) == 1) THEN
  lwx = mx
ELSE
  lwx = 4*mx
END IF
IF (intpol(2) == 1) THEN
  lwy = my+2*mx
ELSE
  lwy = 4*(mx+my)
END IF
IF (lw < lwx+lwy) RETURN
IF (liw < mx+my) RETURN

!     check (xx,yy) grid contained in (x,y) grid

ier = 3
IF (xx(1) < x(1) .OR. xx(mx) > x(nx)) RETURN
IF (yy(1) < y(1) .OR. yy(my) > y(ny)) RETURN

!     check montonicity of grids

ier = 4
DO i=2,nx
  IF (x(i-1) >= x(i)) RETURN
END DO
DO j=2,ny
  IF (y(j-1) >= y(j)) RETURN
END DO
DO ii=2,mx
  IF (xx(ii-1) > xx(ii)) RETURN
END DO
DO jj=2,my
  IF (yy(jj-1) > yy(jj)) RETURN
END DO

!     arguments o.k.

ier = 0

!     set pointer in integer work space

jy = mx+1
IF (intpol(2) == 1) THEN
  
!     linearly interpolate in y
  
  j2 = 1
  j3 = j2
  j4 = j3+my
  j5 = j4
  j6 = j5
  j7 = j6
  j8 = j7+mx
  j9 = j8+mx
  
!     set y interpolation indices and scales and linearly interpolate
  
  CALL linmx(ny,y,my,yy,iw(jy),w(j3))
  i2 = j9
  
!     set work space portion and indices which depend on x interpolation
  
  IF (intpol(1) == 1) THEN
    i3 = i2
    i4 = i3
    i5 = i4
    CALL linmx(nx,x,mx,xx,iw,w(i3))
  ELSE
    i3 = i2+mx
    i4 = i3+mx
    i5 = i4+mx
    CALL cubnmx(nx,x,mx,xx,iw,w(i2),w(i3),w(i4),w(i5))
  END IF
  CALL lint2(nx,ny,p,mx,my,q,intpol,iw(jy),w(j3),  &
      w(j7),w(j8),iw,w(i2),w(i3),w(i4),w(i5))
  RETURN
  
ELSE
  
!     cubically interpolate in y, set indice pointers
  
  j2 = 1
  j3 = j2+my
  j4 = j3+my
  j5 = j4+my
  j6 = j5+my
  j7 = j6+mx
  j8 = j7+mx
  j9 = j8+mx
  CALL cubnmx(ny,y,my,yy,iw(jy),w(j2),w(j3),w(j4),w(j5))
  i2 =  j9+mx
  
!     set work space portion and indices which depend on x interpolation
  
  IF (intpol(1) == 1) THEN
    i3 = i2
    i4 = i3
    i5 = i4
    CALL linmx(nx,x,mx,xx,iw,w(i3))
  ELSE
    i3 = i2+mx
    i4 = i3+mx
    i5 = i4+mx
    CALL cubnmx(nx,x,mx,xx,iw,w(i2),w(i3),w(i4),w(i5))
  END IF
  CALL cubt2(nx,ny,p,mx,my,q,intpol,iw(jy),w(j2),w(j3),  &
      w(j4),w(j5),w(j6),w(j7),w(j8),w(j9),iw,w(i2),w(i3),w(i4),w(i5))
  RETURN
END IF
END SUBROUTINE rgrd2

SUBROUTINE lint2(nx,ny,p,mx,my,q,intpol,jy,dy,pj,pjp, ix,dxm,dx,dxp,dxpp)
IMPLICIT NONE
INTEGER, INTENT(IN )                  :: nx
INTEGER, INTENT(IN )                  :: ny
REAL (KIND=qPrec), INTENT(IN)                     :: p(nx,ny)
INTEGER, INTENT(IN)                      :: mx
INTEGER, INTENT(IN)                      :: my
REAL (KIND=qPrec), INTENT(IN OUT)                        :: q(mx,my)
INTEGER, INTENT(IN)                  :: intpol(2)
INTEGER, INTENT(IN OUT)                      :: jy(my)
REAL (KIND=qPrec), INTENT(IN OUT)                         :: dy(my)
REAL (KIND=qPrec), INTENT(IN OUT)                        :: pj(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                         :: pjp(mx)
INTEGER, INTENT(IN OUT)                  :: ix(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dxm(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dx(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dxp(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dxpp(mx)


INTEGER :: jsave,j,jj,ii




!     linearly interpolate in y

IF (intpol(1) == 1) THEN
  
!     linear in x
  
  jsave = -1
  DO jj=1,my
    j = jy(jj)
    IF (j == jsave) THEN
      
!       j pointer has not moved since last pass (no updates or interpolation)
      
    ELSE IF (j == jsave+1) THEN
      
!       update j and interpolate j+1
      
      DO ii=1,mx
        pj(ii) = pjp(ii)
      END DO
      CALL lint1(nx,p(1,j+1),mx,pjp,ix,dx)
    ELSE
      
!       interpolate j,j+1in pj,pjp on xx mesh
      
      CALL lint1(nx,p(1,j),mx,pj,ix,dx)
      CALL lint1(nx,p(1,j+1),mx,pjp,ix,dx)
    END IF
    
!       save j pointer for next pass
    
    jsave = j
    
!       linearly interpolate q(ii,jj) from pjp,pj in y direction
    
    DO ii=1,mx
      q(ii,jj) = pj(ii)+dy(jj)*(pjp(ii)-pj(ii))
    END DO
  END DO
  
ELSE
  
!     cubic in x
  
  jsave = -1
  DO jj=1,my
    j = jy(jj)
    IF (j == jsave) THEN
      
!       j pointer has not moved since last pass (no updates or interpolation)
      
    ELSE IF (j == jsave+1) THEN
      
!       update j and interpolate j+1
      
      DO ii=1,mx
        pj(ii) = pjp(ii)
      END DO
      CALL cubt1(nx,p(1,j+1),mx,pjp,ix,dxm,dx,dxp,dxpp)
    ELSE
      
!       interpolate j,j+1 in pj,pjp on xx mesh
      
      CALL cubt1(nx,p(1,j),mx,pj,ix,dxm,dx,dxp,dxpp)
      CALL cubt1(nx,p(1,j+1),mx,pjp,ix,dxm,dx,dxp,dxpp)
    END IF
    
!       save j pointer for next pass
    
    jsave = j
    
!       linearly interpolate q(ii,jj) from pjp,pj in y direction
    
    DO ii=1,mx
      q(ii,jj) = pj(ii)+dy(jj)*(pjp(ii)-pj(ii))
    END DO
  END DO
  RETURN
END IF
END SUBROUTINE lint2

SUBROUTINE cubt2(nx,ny,p,mx,my,q,intpol,jy,dym,dy,dyp,  &
    dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp)
IMPLICIT NONE
INTEGER, INTENT(IN )                  :: nx
INTEGER, INTENT(IN )                  :: ny
REAL (KIND=qPrec), INTENT(IN)                     :: p(nx,ny)
INTEGER, INTENT(IN)                      :: mx
INTEGER, INTENT(IN)                      :: my
REAL (KIND=qPrec), INTENT(IN OUT)                        :: q(mx,my)
INTEGER, INTENT(IN)                  :: intpol(2)
INTEGER, INTENT(IN OUT)                      :: jy(my)
REAL (KIND=qPrec), INTENT(IN OUT)                         :: dym(my)
REAL (KIND=qPrec), INTENT(IN OUT)                         :: dy(my)
REAL (KIND=qPrec), INTENT(IN OUT)                         :: dyp(my)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dypp(my)
REAL (KIND=qPrec), INTENT(IN OUT)                        :: pjm(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: pj(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: pjp(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                         :: pjpp(mx)
INTEGER, INTENT(IN OUT)                  :: ix(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dxm(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dx(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dxp(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dxpp(mx)


INTEGER :: jsave,j,jj,ii





IF (intpol(1) == 1) THEN
  
!     linear in x
  
  jsave = -3
  DO jj=1,my
    
!       load closest four j lines containing interpolate on xx mesh
!       for j-1,j,j+1,j+2 in pjm,pj,pjp,pjpp
    
    j = jy(jj)
    IF (j == jsave) THEN
      
!       j pointer has not moved since last pass (no updates or interpolation)
      
    ELSE IF (j == jsave+1) THEN
      
!       update j-1,j,j+1 and interpolate j+2
      
      DO ii=1,mx
        pjm(ii) = pj(ii)
        pj(ii) = pjp(ii)
        pjp(ii) = pjpp(ii)
      END DO
      CALL lint1(nx,p(1,j+2),mx,pjpp,ix,dx)
    ELSE IF (j == jsave+2) THEN
      
!     update j-1,j and interpolate j+1,j+2
      
      DO ii=1,mx
        pjm(ii) = pjp(ii)
        pj(ii) = pjpp(ii)
      END DO
      CALL lint1(nx,p(1,j+1),mx,pjp,ix,dx)
      CALL lint1(nx,p(1,j+2),mx,pjpp,ix,dx)
    ELSE IF (j == jsave+3) THEN
      
!       update j-1 and interpolate j,j+1,j+2
      
      DO ii=1,mx
        pjm(ii) = pjpp(ii)
      END DO
      CALL lint1(nx,p(1,j),mx,pj,ix,dx)
      CALL lint1(nx,p(1,j+1),mx,pjp,ix,dx)
      CALL lint1(nx,p(1,j+2),mx,pjpp,ix,dx)
    ELSE
      
!       interpolate all four j-1,j,j+1,j+2
      
      CALL lint1(nx,p(1,j-1),mx,pjm,ix,dx)
      CALL lint1(nx,p(1,j),mx,pj,ix,dx)
      CALL lint1(nx,p(1,j+1),mx,pjp,ix,dx)
      CALL lint1(nx,p(1,j+2),mx,pjpp,ix,dx)
    END IF
    
!     save j pointer for next pass
    
    jsave = j
    
!     cubically interpolate q(ii,jj) from pjm,pj,pjp,pjpp in y direction
    
    DO ii=1,mx
      q(ii,jj) = dym(jj)*pjm(ii)+dy(jj)*pj(ii)+dyp(jj)*pjp(ii)+  &
          dypp(jj)*pjpp(ii)
    END DO
  END DO
  RETURN
  
ELSE
  
!     cubic in x
  
  jsave = -3
  DO jj=1,my
    
!       load closest four j lines containing interpolate on xx mesh
!       for j-1,j,j+1,j+2 in pjm,pj,pjp,pjpp
    
    j = jy(jj)
    IF (j == jsave) THEN
      
!         j pointer has not moved since last pass (no updates or interpolation)
      
    ELSE IF (j == jsave+1) THEN
      
!         update j-1,j,j+1 and interpolate j+2
      
      DO ii=1,mx
        pjm(ii) = pj(ii)
        pj(ii) = pjp(ii)
        pjp(ii) = pjpp(ii)
      END DO
      CALL cubt1(nx,p(1,j+2),mx,pjpp,ix,dxm,dx,dxp,dxpp)
    ELSE IF (j == jsave+2) THEN
      
!         update j-1,j and interpolate j+1,j+2
      
      DO ii=1,mx
        pjm(ii) = pjp(ii)
        pj(ii) = pjpp(ii)
      END DO
      CALL cubt1(nx,p(1,j+1),mx,pjp,ix,dxm,dx,dxp,dxpp)
      CALL cubt1(nx,p(1,j+2),mx,pjpp,ix,dxm,dx,dxp,dxpp)
    ELSE IF (j == jsave+3) THEN
      
!         update j-1 and interpolate j,j+1,j+2
      
      DO ii=1,mx
        pjm(ii) = pjpp(ii)
      END DO
      CALL cubt1(nx,p(1,j),mx,pj,ix,dxm,dx,dxp,dxpp)
      CALL cubt1(nx,p(1,j+1),mx,pjp,ix,dxm,dx,dxp,dxpp)
      CALL cubt1(nx,p(1,j+2),mx,pjpp,ix,dxm,dx,dxp,dxpp)
    ELSE
      
!         interpolate all four j-1,j,j+1,j+2
      
      CALL cubt1(nx,p(1,j-1),mx,pjm,ix,dxm,dx,dxp,dxpp)
      CALL cubt1(nx,p(1,j),mx,pj,ix,dxm,dx,dxp,dxpp)
      CALL cubt1(nx,p(1,j+1),mx,pjp,ix,dxm,dx,dxp,dxpp)
      CALL cubt1(nx,p(1,j+2),mx,pjpp,ix,dxm,dx,dxp,dxpp)
    END IF
    
!       save j pointer for next pass
    
    jsave = j
    
!       cubically interpolate q(ii,jj) from pjm,pj,pjp,pjpp in y direction
    
    DO ii=1,mx

      q(ii,jj) = dym(jj)*pjm(ii)+dy(jj)*pj(ii)+dyp(jj)*pjp(ii)+  &
          dypp(jj)*pjpp(ii)
    END DO
  END DO
  RETURN
END IF
END SUBROUTINE cubt2

!
 
! Code converted using TO_F90 by Alan Miller
! Date: 2006-11-29  Time: 13:21:12

! ... file rgrd3.f

!     this file contains documentation for subroutine rgrd3 followed by
!     fortran code for rgrd3 and additional subroutines.

! ... author

!     John C. Adams (NCAR 1997)

! ... subroutine rgrd3(nx,ny,nz,x,y,z,p,mx,my,mz,xx,yy,zz,q,intpol,
!    +                 w,lw,iw,liw,ier)

! ... purpose

!     subroutine rgrd3 interpolates the values p(i,j,k) on the orthogonal
!     grid (x(i),y(j),z(k)) for i=1,...,nx; j=1,...,ny; k=1,...,nz
!     onto q(ii,jj,kk) on the orthogonal grid (xx(ii),yy(jj),zz(kk)) for
!     ii=1,...,mx; jj=1,...,my; kk=1,...,mz.

! ... language

!     coded in portable FORTRAN77 and FORTRAN90

! ... test program

!     file trgrd3.f on regridpack includes a test program for subroutine rgrd3

! ... method

!     linear or cubic interpolation is used (independently) in
!     each direction (see argument intpol).

! ... required files

!     files rgrd2.f and rgrd1.f must be loaded with rgrd3.f.  they
!     include subroutines called by the routines in rgrd3.f

! ... requirements

!     each of the x,y,z grids must be strictly montonically increasing
!     and each of the xx,yy,zz grids must be montonically increasing
!     (see ier = 4).  in addition the (X,Y,Z) region

!          [xx(1),xx(mx)] X [yy(1),yy(my)] X [zz(1),zz(mz)]

!     must lie within the (X,Y,Z) region

!          [x(1),x(nx)] X [y(1),y(ny)] X [z(1),z(nz)].

!     extrapolation is not allowed (see ier=3).  if these (X,Y,Z)
!     regions are identical and the orthogonal grids are UNIFORM
!     in each direction then subroutine rgrd3u (see file rgrd3u.f)
!     should be used instead of rgrd3.

! ... efficiency

!     inner most loops in regridpack software vectorize.  if
!     the integer arguments mx,my,mz (see below) have different values,
!     optimal vectorization will be achieved if mx > my > mz.

! *** input arguments

! ... nx

!     the integer dimension of the grid vector x and the first dimension of p.
!     nx > 1 if intpol(1) = 1 or nx > 3 if intpol(1) = 3 is required.

! ... ny

!     the integer dimension of the grid vector y and the second dimension of p.
!     ny > 1 if intpol(2) = 1 or ny > 3 if intpol(2) = 3 is required.

! ... nz

!     the integer dimension of the grid vector z and the third dimension of p.
!     nz > 1 if intpol(3) = 1 or nz > 3 if intpol(3) = 3 is required.

! ... x

!     a real nx vector of strictly increasing values which defines the x
!     portion of the orthogonal grid on which p is given

! ... y

!     a real ny vector of strictly increasing values which defines the y
!     portion of the orthogonal grid on which p is given

! ... z

!     a real nz vector of strictly increasing values which defines the z
!     portion of the orthogonal grid on which p is given

! ... p

!     a real nx by ny by nz array of values given on the (x,y,z) grid

! ... mx

!     the integer dimension of the grid vector xx and the first dimension of q.
!     mx > 0 is required.

! ... my

!     the integer dimension of the grid vector yy and the second dimension of q.
!     my > 0 is required.

! ... mz

!     the integer dimension of the grid vector zz and the third dimension of q.
!     mz > 0 is required.

! ... xx

!     a real mx vector of increasing values which defines the x portion of the
!     orthogonal grid on which q is defined.  xx(1) < x(1) or xx(mx) > x(nx)
!     is not allowed (see ier = 3)

! ... yy

!     a real my vector of increasing values which defines the y portion of the
!     orthogonal grid on which q is defined.  yy(1) < y(1) or yy(my) > y(ny)
!     is not allowed (see ier = 3)

! ... zz

!     a real mz vector of increasing values which defines the z portion of the
!     orthogonal grid on which q is defined.  zz(1) < z(1) or zz(mz) > z(nz)
!     is not allowed (see ier = 3)

! ... intpol

!     an integer vector of dimension 3 which sets linear or cubic
!     interpolation in each of the x,y,z directions as follows:

!        intpol(1) = 1 sets linear interpolation in the x direction
!        intpol(1) = 3 sets cubic interpolation in the x direction.

!        intpol(2) = 1 sets linear interpolation in the y direction
!        intpol(2) = 3 sets cubic interpolation in the y direction.

!        intpol(3) = 1 sets linear interpolation in the z direction
!        intpol(3) = 3 sets cubic interpolation in the z direction.

!     values other than 1 or 3 in intpol are not allowed (ier = 5).

! ... w

!     a real work space of length at least lw which must be provided in the
!     routine calling rgrd3


! ... lw

!     the integer length of the real work space w.  let

!          lwx =   mx            if intpol(1) = 1
!          lwx = 4*mx            if intpol(1) = 3

!          lwy = my+2*mx         if intpol(2) = 1
!          lwy = 4*(mx+my)       if intpol(2) = 3

!          lwz = 2*mx*my+mz      if intpol(3) = 1
!          lwz = 4*(mx*my+mz)    if intpol(3) = 3

!     then lw must be greater than or equal to lwx+lwy+lwz

! ... iw

!     an integer work space of length at least liw which must be provided in the
!     routine calling rgrd3

! ... liw

!     the integer length of the integer work space iw.  liw must be at least mx+my+mz

! *** output arguments


! ... q

!     a real mx by my by mz array of values on the (xx,yy,zz) grid which are
!     interpolated from p on the (x,y,z) grid

! ... ier

!     an integer error flag set as follows:

!     ier = 0 if no errors in input arguments are detected

!     ier = 1 if  min0(mx,my,mz) < 1

!     ier = 2 if nx < 2 when intpol(1)=1 or nx < 4 when intpol(1)=3 (or)
!                ny < 2 when intpol(2)=1 or ny < 4 when intpol(2)=3 (or)
!                nz < 2 when intpol(3)=1 or nz < 4 when intpol(3)=3.

!     ier = 3 if xx(1) < x(1) or x(nx) < xx(mx) (or)
!                yy(1) < y(1) or y(ny) < yy(my) (or)
!                zz(1) < z(1) or z(nz) < zz(mz)

! *** to avoid this flag when end points are intended to be the
!     same but may differ slightly due to roundoff error, they
!     should be set exactly in the calling routine (e.g., if both
!     grids have the same y boundaries then yy(1)=y(1) and yy(my)=y(ny)
!     should be set before calling rgrd3)

!     ier = 4 if one of the grids x,y,z is not strictly monotonically
!             increasing or if one of the grids xx,yy,zz is not
!             montonically increasing.  more precisely if:

!             x(i+1) <= x(i) for some i such that 1 <= i < nx (or)

!             y(j+1) <= y(j) for some j such that 1 <= j < ny (or)

!             z(k+1) <= z(k) for some k such that 1 <= k < nz (or)

!             xx(ii+1) < xx(ii) for some ii such that 1 <= ii < mx (or)

!             yy(jj+1) < yy(jj) for some jj such that 1 <= jj < my (or)

!             zz(kk+1) < zz(k)  for some kk such that 1 <= kk < mz

!     ier = 5 if lw or liw is too small (insufficient work space)

!     ier = 6 if any of intpol(1),intpol(2),intpol(3) is not equal to 1 or 3

! ************************************************************************

!     end of rgrd3 documentation, fortran code follows:

! ************************************************************************

SUBROUTINE rgrd3(nx,ny,nz,x,y,z,p,mx,my,mz,xx,yy,zz,q,intpol, w,lw,iw,liw,ier)
IMPLICIT NONE
INTEGER, INTENT(IN)                      :: nx
INTEGER, INTENT(IN)                      :: ny
INTEGER, INTENT(IN)                      :: nz
REAL (KIND=qPrec), INTENT(IN)                         :: x(nx)
REAL (KIND=qPrec), INTENT(IN)                         :: y(ny)
REAL (KIND=qPrec), INTENT(IN)                         :: z(nz)
REAL (KIND=qPrec), INTENT(IN)                     :: p(nx,ny,nz)
INTEGER, INTENT(IN)                      :: mx
INTEGER, INTENT(IN)                      :: my
INTEGER, INTENT(IN)                      :: mz
REAL (KIND=qPrec), INTENT(IN)                         :: xx(mx)
REAL (KIND=qPrec), INTENT(IN)                         :: yy(my)
REAL (KIND=qPrec), INTENT(IN)                         :: zz(mz)
REAL (KIND=qPrec), INTENT(OUT)                     :: q(mx,my,mz)
INTEGER, INTENT(IN)                      :: intpol(3)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: w(lw)
INTEGER, INTENT(IN OUT)                  :: lw
INTEGER, INTENT(IN OUT)                  :: iw(liw)
INTEGER, INTENT(IN OUT)                  :: liw
INTEGER, INTENT(IN OUT)                     :: ier


INTEGER :: lwx,lwy,lwz,jy,kz,mxmy




INTEGER :: i,ii,j,jj,k,kk
INTEGER :: i2,i3,i4,i5
INTEGER :: j2,j3,j4,j5,j6,j7,j8,j9
INTEGER :: k2,k3,k4,k5,k6,k7,k8,k9

!     check input arguments

ier = 1

!     check (xx,yy,zz) grid resolution

IF (MIN0(mx,my,mz) < 1) RETURN

!     check intpol

ier = 6
IF (intpol(1) /= 1 .AND. intpol(1) /= 3) RETURN
IF (intpol(2) /= 1 .AND. intpol(2) /= 3) RETURN
IF (intpol(3) /= 1 .AND. intpol(3) /= 3) RETURN

!     check (x,y,z) grid resolution

ier = 2
IF (intpol(1) == 1 .AND. nx < 2) RETURN
IF (intpol(1) == 3 .AND. nx < 4) RETURN
IF (intpol(2) == 1 .AND. ny < 2) RETURN
IF (intpol(2) == 3 .AND. ny < 4) RETURN
IF (intpol(3) == 1 .AND. nz < 2) RETURN
IF (intpol(3) == 3 .AND. nz < 4) RETURN

!     check work space length input and set minimum

ier = 5
mxmy = mx*my
IF (intpol(1) == 1) THEN
  lwx = mx
ELSE
  lwx = 4*mx
END IF
IF (intpol(2) == 1) THEN
  lwy = (my+2*mx)
ELSE
  lwy = 4*(my+mx)
END IF
IF (intpol(3) == 1) THEN
  lwz = (2*mxmy+mz)
ELSE
  lwz = 4*(mxmy+mz)
END IF
IF (lw < lwx+lwy+lwz) RETURN
IF (liw < mx+my+mz) RETURN

!     check (xx,yy,zz) grid contained in (x,y,z) grid

ier = 3
IF (xx(1) < x(1) .OR. xx(mx) > x(nx)) RETURN
IF (yy(1) < y(1) .OR. yy(my) > y(ny)) RETURN
IF (zz(1) < z(1) .OR. zz(mz) > z(nz)) RETURN

!     check montonicity of grids

ier = 4
DO i=2,nx
  IF (x(i-1) >= x(i)) RETURN
END DO
DO j=2,ny
  IF (y(j-1) >= y(j)) RETURN
END DO
DO k=2,nz
  IF (z(k-1) >= z(k)) RETURN
END DO
DO ii=2,mx
  IF (xx(ii-1) > xx(ii)) RETURN
END DO
DO jj=2,my
  IF (yy(jj-1) > yy(jj)) RETURN
END DO
DO kk=2,mz
  IF (zz(kk-1) > zz(kk)) RETURN
END DO

!     arguments o.k.

ier = 0
jy = mx+1
kz = mx+my+1
IF (intpol(3) == 1) THEN
  
!     linearly interpolate in nz, set work space pointers and scales
  
  k2 = 1
  k3 = k2
  k4 = k3+mz
  k5 = k4
  k6 = k5
  k7 = k6
  k8 = k7+mxmy
  k9 = k8+mxmy
  CALL linmx(nz,z,mz,zz,iw(kz),w(k3))
  j2 = k9
  
!     set indices and scales which depend on y interpolation
  
  IF (intpol(2) == 1) THEN
!     linear in y
    j3 = j2
    j4 = j3+my
    j5 = j4
    j6 = j5
    j7 = j6
    j8 = j7+mx
    j9 = j8+mx
    CALL linmx(ny,y,my,yy,iw(jy),w(j3))
    i2 = j9
  ELSE
!     cubic in y
    j3 = j2+my
    j4 = j3+my
    j5 = j4+my
    j6 = j5+my
    j7 = j6+mx
    j8 = j7+mx
    j9 = j8+mx
    CALL cubnmx(ny,y,my,yy,iw(jy),w(j2),w(j3),w(j4),w(j5))
    i2 = j9+mx
  END IF
  
!     set indices and scales which depend on x interpolation
  
  IF (intpol(1) == 1) THEN
!     linear in x
    i3 = i2
    i4 = i3
    i5 = i4
    CALL linmx(nx,x,mx,xx,iw,w(i3))
  ELSE
!     cubic in x
    i3 = i2+mx
    i4 = i3+mx
    i5 = i4+mx
    CALL cubnmx(nx,x,mx,xx,iw,w(i2),w(i3),w(i4),w(i5))
  END IF
  CALL lint3(nx,ny,nz,p,mx,my,mxmy,mz,q,intpol,iw(kz),  &
      w(k3),w(k7),w(k8),iw(jy),w(j2),w(j3),w(j4),w(j5),w(j6),  &
      w(j7),w(j8),w(j9),iw,w(i2),w(i3),w(i4),w(i5))
  RETURN
ELSE
  
!     cubically interpolate in z
  
  k2 = 1
  k3 = k2+mz
  k4 = k3+mz
  k5 = k4+mz
  k6 = k5+mz
  k7 = k6+mxmy
  k8 = k7+mxmy
  k9 = k8+mxmy
  CALL cubnmx(nz,z,mz,zz,iw(kz),w(k2),w(k3),w(k4),w(k5))
  j2 = k9+mxmy
  
!     set indices which depend on y interpolation
  
  IF (intpol(2) == 1) THEN
    j3 = j2
    j4 = j3+my
    j5 = j4
    j6 = j5
    j7 = j6
    j8 = j7+mx
    j9 = j8+mx
    CALL linmx(ny,y,my,yy,iw(jy),w(j3))
    i2 = j9
  ELSE
    j3 = j2+my
    j4 = j3+my
    j5 = j4+my
    j6 = j5+my
    j7 = j6+mx
    j8 = j7+mx
    j9 = j8+mx
    CALL cubnmx(ny,y,my,yy,iw(jy),w(j2),w(j3),w(j4),w(j5))
    i2 = j9+mx
  END IF
  
!     set work space portion and indices which depend on x interpolation
  
  IF (intpol(1) == 1) THEN
    i3 = i2
    i4 = i3
    i5 = i4
    CALL linmx(nx,x,mx,xx,iw,w(i3))
  ELSE
    i3 = i2+mx
    i4 = i3+mx
    i5 = i4+mx
    CALL cubnmx(nx,x,mx,xx,iw,w(i2),w(i3),w(i4),w(i5))
  END IF
  CALL cubt3(nx,ny,nz,p,mx,my,mxmy,mz,q,intpol,  &
      iw(kz),w(k2),w(k3),w(k4),w(k5),w(k6),w(k7),w(k8),w(k9),  &
      iw(jy),w(j2),w(j3),w(j4),w(j5),w(j6),w(j7),w(j8),w(j9),  &
      iw,w(i2),w(i3),w(i4),w(i5))
  RETURN
  
END IF

END SUBROUTINE rgrd3

SUBROUTINE lint3(nx,ny,nz,p,mx,my,mxmy,mz,q,intpol,kz,  &
    dz,pk,pkp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp)
IMPLICIT NONE
!     linearly interpolate in z direction


INTEGER, INTENT(IN)                  :: nx
INTEGER, INTENT(IN)                  :: ny
INTEGER, INTENT(IN)                  :: nz
REAL (KIND=qPrec), INTENT(IN)                     :: p(nx,ny,nz)
INTEGER, INTENT(IN)                  :: mx
INTEGER, INTENT(IN)                  :: my
INTEGER, INTENT(IN)                      :: mxmy
INTEGER, INTENT(IN)                      :: mz
REAL (KIND=qPrec), INTENT(OUT)                        :: q(mxmy,mz)
INTEGER, INTENT(IN)                  :: intpol(3)
INTEGER, INTENT(IN OUT)                      :: kz(mz)
REAL (KIND=qPrec), INTENT(IN OUT)                         :: dz(mz)
REAL (KIND=qPrec), INTENT(IN OUT)                        :: pk(mxmy)
REAL (KIND=qPrec), INTENT(IN OUT)                         :: pkp(mxmy)
INTEGER, INTENT(IN OUT)                  :: jy(my)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dym(my)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dy(my)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dyp(my)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dypp(my)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: pjm(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: pj(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: pjp(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: pjpp(mx)
INTEGER, INTENT(IN OUT)                  :: ix(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dxm(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dx(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dxp(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dxpp(mx)








INTEGER :: k,kk,iijj,ksave

IF (intpol(2) == 1) THEN
  
!     linear in y
  
  ksave = -1
  DO kk=1,mz
    k = kz(kk)
    IF (k == ksave) THEN
      
!     k pointer has not moved since last pass (no updates or interpolation)
      
    ELSE IF (k == ksave+1) THEN
      
!     update k and interpolate k+1
      
      DO iijj=1,mxmy
        pk(iijj) = pkp(iijj)
      END DO
      CALL lint2(nx,ny,p(1,1,k+1),mx,my,pkp,intpol,jy,dy,  &
          pj,pjp,ix,dxm,dx,dxp,dxpp)
    ELSE
      
!     interpolate k,k+1 in pk,pkp on xx,yy mesh
      
      CALL lint2(nx,ny,p(1,1,k),mx,my,pk,intpol,jy,dy,  &
          pj,pjp,ix,dxm,dx,dxp,dxpp)
      CALL lint2(nx,ny,p(1,1,k+1),mx,my,pkp,intpol,jy,dy,  &
          pj,pjp,ix,dxm,dx,dxp,dxpp)
    END IF
    
!     save k pointer for next pass
    
    ksave = k
    
!     linearly interpolate q(ii,jj,k) from pk,pkp in z direction
    
    DO iijj=1,mxmy
      q(iijj,kk) = pk(iijj)+dz(kk)*(pkp(iijj)-pk(iijj))
    END DO
  END DO
  RETURN
  
ELSE
  
!     cubic in y
  
  ksave = -1
  DO kk=1,mz
    k = kz(kk)
    IF (k == ksave) THEN
      
!     k pointer has not moved since last pass (no updates or interpolation)
      
    ELSE IF (k == ksave+1) THEN
      
!     update k and interpolate k+1
      
      DO iijj=1,mxmy
        pk(iijj) = pkp(iijj)
      END DO
      CALL cubt2(nx,ny,p(1,1,k+1),mx,my,pkp,intpol,  &
          jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp)
    ELSE
      
!     interpolate k,k+1 in pk,pkp on xx,yy mesh
      
      CALL cubt2(nx,ny,p(1,1,k),mx,my,pk,intpol,  &
          jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp)
      CALL cubt2(nx,ny,p(1,1,k+1),mx,my,pkp,intpol,  &
          jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp)
    END IF
    
!     save k pointer for next pass
    
    ksave = k
    
!     linearly interpolate q(ii,jj,k) from pk,pkp in z direction
    
    DO iijj=1,mxmy
      q(iijj,kk) = pk(iijj)+dz(kk)*(pkp(iijj)-pk(iijj))
    END DO
  END DO
  RETURN
  
END IF
END SUBROUTINE lint3

SUBROUTINE cubt3(nx,ny,nz,p,mx,my,mxmy,mz,q,intpol,  &
    kz,dzm,dz,dzp,dzpp,pkm,pk,pkp,pkpp,jy,dym,dy,dyp,dypp,pjm,pj,  &
    pjp,pjpp,ix,dxm,dx,dxp,dxpp)

!     cubically interpolate in z

IMPLICIT NONE
INTEGER, INTENT(IN)                  :: nx
INTEGER, INTENT(IN)                  :: ny
INTEGER, INTENT(IN)                  :: nz
REAL (KIND=qPrec), INTENT(IN)                     :: p(nx,ny,nz)
INTEGER, INTENT(IN)                  :: mx
INTEGER, INTENT(IN)                  :: my
INTEGER, INTENT(IN OUT)                      :: mxmy
INTEGER, INTENT(IN)                      :: mz
REAL (KIND=qPrec), INTENT(IN OUT)                        :: q(mxmy,mz)
INTEGER, INTENT(IN)                  :: intpol(3)
INTEGER, INTENT(IN OUT)                      :: kz(mz)
REAL (KIND=qPrec), INTENT(IN OUT)                         :: dzm(mz)
REAL (KIND=qPrec), INTENT(IN OUT)                         :: dz(mz)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dzp(mz)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dzpp(mz)
REAL (KIND=qPrec), INTENT(IN OUT)                        :: pkm(mxmy)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: pk(mxmy)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: pkp(mxmy)
REAL (KIND=qPrec), INTENT(IN OUT)                         :: pkpp(mxmy)
INTEGER, INTENT(IN OUT)                  :: jy(my)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dym(my)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dy(my)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dyp(my)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dypp(my)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: pjm(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: pj(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: pjp(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: pjpp(mx)
INTEGER, INTENT(IN OUT)                  :: ix(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dxm(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dx(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dxp(mx)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: dxpp(mx)

INTEGER :: k,kk,ksave,iijj








IF (intpol(2) == 1) THEN
  
!       linear in y
  
  ksave = -3
  DO kk=1,mz
    k = kz(kk)
    IF (k == ksave) THEN
      
!       k pointer has not moved since last pass (no updates or interpolation)
      
    ELSE IF (k == ksave+1) THEN
      
!       update k-1,k,k+1 and interpolate k+2
      
      DO iijj=1,mxmy
        pkm(iijj) = pk(iijj)
        pk(iijj) = pkp(iijj)
        pkp(iijj) = pkpp(iijj)
      END DO
      CALL lint2(nx,ny,p(1,1,k+2),mx,my,pkpp,  &
          intpol,jy,dy,pj,pjp,ix,dxm,dx,dxp,dxpp)
    ELSE IF (k == ksave+2) THEN
      
!       update k-1,k and interpolate k+1,k+2
      
      DO iijj=1,mxmy
        pkm(iijj) = pkp(iijj)
        pk(iijj) = pkpp(iijj)
      END DO
      CALL lint2(nx,ny,p(1,1,k+1),mx,my,pkp,  &
          intpol,jy,dy,pj,pjp,ix,dxm,dx,dxp,dxpp)
      CALL lint2(nx,ny,p(1,1,k+2),mx,my,pkpp,  &
          intpol,jy,dy,pj,pjp,ix,dxm,dx,dxp,dxpp)
    ELSE IF (k == ksave+3) THEN
      
!       update k-1 and interpolate k,k+1,k+2
      
      DO iijj=1,mxmy
        pkm(iijj) = pkpp(iijj)
      END DO
      CALL lint2(nx,ny,p(1,1,k),mx,my,pk,  &
          intpol,jy,dy,pj,pjp,ix,dxm,dx,dxp,dxpp)
      CALL lint2(nx,ny,p(1,1,k+1),mx,my,pkp,  &
          intpol,jy,dy,pj,pjp,ix,dxm,dx,dxp,dxpp)
      CALL lint2(nx,ny,p(1,1,k+2),mx,my,pkpp,  &
          intpol,jy,dy,pj,pjp,ix,dxm,dx,dxp,dxpp)
    ELSE
      
!       interpolate all four k-1,k,k+1,k+2
      
      CALL lint2(nx,ny,p(1,1,k-1),mx,my,pkm,  &
          intpol,jy,dy,pj,pjp,ix,dxm,dx,dxp,dxpp)
      CALL lint2(nx,ny,p(1,1,k),mx,my,pk,  &
          intpol,jy,dy,pj,pjp,ix,dxm,dx,dxp,dxpp)
      CALL lint2(nx,ny,p(1,1,k+1),mx,my,pkp,  &
          intpol,jy,dy,pj,pjp,ix,dxm,dx,dxp,dxpp)
      CALL lint2(nx,ny,p(1,1,k+2),mx,my,pkpp,  &
          intpol,jy,dy,pj,pjp,ix,dxm,dx,dxp,dxpp)
    END IF
    
!       save k pointer for next pass
    
    ksave = k
    
!       cubically interpolate q(ii,jj,kk) from pkm,pk,pkp,pkpp in z direction
    
    DO iijj=1,mxmy
      q(iijj,kk) = dzm(kk)*pkm(iijj) + dz(kk)*pk(iijj) +  &
          dzp(kk)*pkp(iijj) + dzpp(kk)*pkpp(iijj)
    END DO
  END DO
  RETURN
  
ELSE
  
!       cubic in y
  
  ksave = -3
  DO kk=1,mz
    k = kz(kk)
    IF (k == ksave) THEN
      
!       k pointer has not moved since last pass (no updates or interpolation)
      
    ELSE IF (k == ksave+1) THEN
      
!       update k-1,k,k+1 and interpolate k+2
      
      DO iijj=1,mxmy
        pkm(iijj) = pk(iijj)
        pk(iijj) = pkp(iijj)
        pkp(iijj) = pkpp(iijj)
      END DO
      CALL cubt2(nx,ny,p(1,1,k+2),mx,my,pkpp,intpol,jy,dym,dy,  &
          dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp)
    ELSE IF (k == ksave+2) THEN
      
!       update k-1,k and interpolate k+1,k+2
      
      DO iijj=1,mxmy
        pkm(iijj) = pkp(iijj)
        pk(iijj) = pkpp(iijj)
      END DO
      CALL cubt2(nx,ny,p(1,1,k+1),mx,my,pkp,intpol,jy,dym,dy,  &
          dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp)
      CALL cubt2(nx,ny,p(1,1,k+2),mx,my,pkpp,intpol,jy,dym,dy,  &
          dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp)
    ELSE IF (k == ksave+3) THEN
      
!       update k-1 and interpolate k,k+1,k+2
      
      DO iijj=1,mxmy
        pkm(iijj) = pkpp(iijj)
      END DO
      CALL cubt2(nx,ny,p(1,1,k),mx,my,pk,intpol,  &
          jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp)
      CALL cubt2(nx,ny,p(1,1,k+1),mx,my,pkp,intpol,  &
          jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp)
      CALL cubt2(nx,ny,p(1,1,k+2),mx,my,pkpp,intpol,  &
          jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp)
    ELSE
      
!     interpolate all four k-1,k,k+1,k+2
      
      CALL cubt2(nx,ny,p(1,1,k-1),mx,my,pkm,intpol,  &
          jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp)
      CALL cubt2(nx,ny,p(1,1,k),mx,my,pk,intpol,  &
          jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp)
      CALL cubt2(nx,ny,p(1,1,k+1),mx,my,pkp,intpol,  &
          jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp)
      CALL cubt2(nx,ny,p(1,1,k+2),mx,my,pkpp,intpol,  &
          jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp)
      
    END IF
    
!       save k pointer for next pass
    
    ksave = k
    
!       cubically interpolate q(ii,jj,kk) from pkm,pk,pkp,pkpp in z direction
    
    DO iijj=1,mxmy
      q(iijj,kk) = dzm(kk)*pkm(iijj) + dz(kk)*pk(iijj) +  &
          dzp(kk)*pkp(iijj) + dzpp(kk)*pkpp(iijj)
    END DO
  END DO
  RETURN
END IF

END SUBROUTINE cubt3

!
 
! Code converted using TO_F90 by Alan Miller
! Date: 2006-11-29  Time: 13:20:43
 
!     file mudcom.f

!  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
!  .                                                             .
!  .                  copyright (c) 1999 by UCAR                 .
!  .                                                             .
!  .       UNIVERSITY CORPORATION for ATMOSPHERIC RESEARCH       .
!  .                                                             .
!  .                      all rights reserved                    .
!  .                                                             .
!  .                                                             .
!  .                      MUDPACK version 5.0                    .
!  .                                                             .
!  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

! ... author and specialist

!          John C. Adams (National Center for Atmospheric Research)
!          email: johnad@ucar.edu, phone: 303-497-1213

! ... For MUDPACK 5.0 information, visit the website:
!     (http://www.scd.ucar.edu/css/software/mudpack)

! ... purpose

!     mudcom.f is a common subroutines file containing subroutines
!     called by some or all of the real two- and three-dimensional
!     mudpack solvers.  mudcom.f must be loaded with any real mudpack
!     solver.

SUBROUTINE swk2(nfx,nfy,phif,rhsf,phi,rhs)

!     set phif,rhsf input in arrays which include
!     virtual boundaries for phi (for all 2-d real codes)

IMPLICIT NONE
INTEGER, INTENT(IN OUT)                      :: nfx
INTEGER, INTENT(IN OUT)                  :: nfy
REAL (KIND=qPrec), INTENT(IN OUT)                         :: phif(nfx,nfy)
REAL (KIND=qPrec), INTENT(IN OUT)                         :: rhsf(nfx,nfy)
REAL (KIND=qPrec), INTENT(IN OUT)                        :: phi(0:nfx+1,0:nfy+1)
REAL (KIND=qPrec), INTENT(IN OUT)                        :: rhs(nfx,nfy)

INTEGER :: i,j



DO j=1,nfy
  DO i=1,nfx
    phi(i,j) = phif(i,j)
    rhs(i,j) = rhsf(i,j)
  END DO
END DO

!     set virtual boundaries in phi to zero

DO j=0,nfy+1
  phi(0,j) = 0.0
  phi(nfx+1,j) = 0.0
END DO
DO i=0,nfx+1
  phi(i,0) = 0.0
  phi(i,nfy+1) = 0.0
END DO
RETURN
END SUBROUTINE swk2

SUBROUTINE trsfc2(nx,ny,phi,rhs,ncx,ncy,phic,rhsc)

!     transfer fine grid to coarse grid

IMPLICIT NONE
INTEGER, INTENT(IN OUT)                      :: nx
INTEGER, INTENT(IN OUT)                  :: ny
REAL (KIND=qPrec), INTENT(IN OUT)                         :: phi(0:nx+1,0:ny+1)
REAL (KIND=qPrec), INTENT(IN OUT)                         :: rhs(nx,ny)
INTEGER, INTENT(IN OUT)                      :: ncx
INTEGER, INTENT(IN OUT)                  :: ncy
REAL (KIND=qPrec), INTENT(IN OUT)                        :: phic(0:ncx+1,0:ncy+1)
REAL (KIND=qPrec), INTENT(IN OUT)                        :: rhsc(ncx,ncy)

INTEGER :: i,j,ic,jc



!     set virtual boundaries in phic to zero

DO jc=0,ncy+1
  phic(0,jc) = 0.0
  phic(ncx+1,jc) = 0.0
END DO
DO ic=0,ncx+1
  phic(ic,0) = 0.0
  phic(ic,ncy+1) = 0.0
END DO
IF (ncx < nx .AND. ncy < ny) THEN
  
!     coarsening in both x and y
  
  DO jc=1,ncy
    j = jc+jc-1
    DO ic=1,ncx
      i = ic+ic-1
      phic(ic,jc) = phi(i,j)
      rhsc(ic,jc) = rhs(i,j)
    END DO
  END DO
ELSE IF (ncx < nx .AND. ncy == ny) THEN
  
!     coarsening in x only
  
  DO jc=1,ncy
    j = jc
    DO ic=1,ncx
      i = ic+ic-1
      phic(ic,jc) = phi(i,j)
      rhsc(ic,jc) = rhs(i,j)
    END DO
  END DO
ELSE
  
!     coarsening in y only
  
  DO jc=1,ncy
    j = jc+jc-1
    DO ic=1,ncx
      i = ic
      phic(ic,jc) = phi(i,j)
      rhsc(ic,jc) = rhs(i,j)
    END DO
  END DO
END IF
RETURN
END SUBROUTINE trsfc2

SUBROUTINE res2(nx,ny,resf,ncx,ncy,rhsc,nxa,nxb,nyc,nyd)
IMPLICIT NONE
INTEGER, INTENT(IN)                      :: nx
INTEGER, INTENT(IN)                      :: ny
REAL (KIND=qPrec), INTENT(IN OUT)                         :: resf(nx,ny)
INTEGER, INTENT(IN )                      :: ncx
INTEGER, INTENT(IN )                      :: ncy
REAL (KIND=qPrec), INTENT(IN OUT)                        :: rhsc(ncx,ncy)
INTEGER, INTENT(IN OUT)                      :: nxa
INTEGER, INTENT(IN OUT)                  :: nxb
INTEGER, INTENT(IN OUT)                      :: nyc
INTEGER, INTENT(IN OUT)                  :: nyd


INTEGER :: i,j,ic,jc,im1,ip1,jm1,jp1,ix,jy

!     restrict fine grid residual in resf to coarse grid in rhsc
!     using full weighting for all real 2d codes



!     set x,y coarsening integer subscript scales

ix = 1
IF (ncx == nx) ix = 0
jy = 1
IF (ncy == ny) jy = 0

!     restrict on interior

IF (ncy < ny .AND. ncx < nx) THEN
  
!     coarsening in both directions
  
!$OMP PARALLEL DO PRIVATE(i,j,ic,jc), SHARED(resf,rhsc,ncx,ncy)
  DO jc=2,ncy-1
    j = jc+jc-1
    DO ic=2,ncx-1
      i = ic+ic-1
      rhsc(ic,jc) = (resf(i-1,j-1)+resf(i+1,j-1)+resf(i-1,j+1)+  &
          resf(i+1,j+1)+2.*(resf(i-1,j)+resf(i+1,j)+  &
          resf(i,j-1)+resf(i,j+1))+4.*resf(i,j))*.0625
    END DO
  END DO
ELSE IF (ncy == ny) THEN
  
!     no coarsening in y but coarsening in x
  
!$OMP PARALLEL DO PRIVATE(i,j,ic,jc), SHARED(resf,rhsc,ncx,ncy)
  DO jc=2,ncy-1
    j = jc
    DO ic=2,ncx-1
      i = ic+ic-1
      rhsc(ic,jc) = (resf(i-1,j-1)+resf(i+1,j-1)+resf(i-1,j+1)+  &
          resf(i+1,j+1)+2.*(resf(i-1,j)+resf(i+1,j)+  &
          resf(i,j-1)+resf(i,j+1))+4.*resf(i,j))*.0625
    END DO
  END DO
ELSE
  
!     no coarsening in x but coarsening in y
  
!$OMP PARALLEL DO PRIVATE(i,j,ic,jc), SHARED(resf,rhsc,ncx,ncy)
  DO jc=2,ncy-1
    j = jc+jc-1
    DO ic=2,ncx-1
      i = ic
      rhsc(ic,jc) = (resf(i-1,j-1)+resf(i+1,j-1)+resf(i-1,j+1)+  &
          resf(i+1,j+1)+2.*(resf(i-1,j)+resf(i+1,j)+  &
          resf(i,j-1)+resf(i,j+1))+4.*resf(i,j))*.0625
    END DO
  END DO
END IF

!     set residual on boundaries

DO jc=1,ncy,ncy-1
  
!     y=yc,yd boundaries
  
  j = jc+jy*(jc-1)
  jm1 = MAX0(j-1,2)
  jp1 = MIN0(j+1,ny-1)
  IF (j == 1 .AND. nyc == 0) jm1 = ny-1
  IF (j == ny .AND. nyc == 0) jp1 = 2
  
!     y=yc,yd and x=xa,xb cornors
  
  DO ic=1,ncx,ncx-1
    i = ic+ix*(ic-1)
    im1 = MAX0(i-1,2)
    ip1 = MIN0(i+1,nx-1)
    IF (i == 1 .AND. nxa == 0) im1 = nx-1
    IF (i == nx .AND. nxa == 0) ip1 = 2
    rhsc(ic,jc) = (resf(im1,jm1)+resf(ip1,jm1)+resf(im1,jp1)+  &
        resf(ip1,jp1)+2.*(resf(im1,j)+resf(ip1,j)+  &
        resf(i,jm1)+resf(i,jp1))+4.*resf(i,j))*.0625
  END DO
  
!     set y=yc,yd interior edges
  
  DO ic=2,ncx-1
    i = ic+ix*(ic-1)
    rhsc(ic,jc) = (resf(i-1,jm1)+resf(i+1,jm1)+resf(i-1,jp1)+  &
        resf(i+1,jp1)+2.*(resf(i-1,j)+resf(i+1,j)+  &
        resf(i,jm1)+resf(i,jp1))+4.*resf(i,j))*.0625
  END DO
END DO

!     set x=xa,xb interior edges

DO ic=1,ncx,ncx-1
  i = ic+ix*(ic-1)
  im1 = MAX0(i-1,2)
  ip1 = MIN0(i+1,nx-1)
  IF (i == 1 .AND. nxa == 0) im1 = nx-1
  IF (i == nx .AND. nxa == 0) ip1 = 2
  DO jc=2,ncy-1
    j = jc+jy*(jc-1)
    rhsc(ic,jc) = (resf(im1,j-1)+resf(ip1,j-1)+resf(im1,j+1)+  &
        resf(ip1,j+1)+2.*(resf(im1,j)+resf(ip1,j)+  &
        resf(i,j-1)+resf(i,j+1))+4.*resf(i,j))*.0625
  END DO
END DO

!     set coarse grid residual zero on specified boundaries

IF (nxa == 1) THEN
  DO jc=1,ncy
    rhsc(1,jc) = 0.0
  END DO
END IF
IF (nxb == 1) THEN
  DO jc=1,ncy
    rhsc(ncx,jc) = 0.0
  END DO
END IF
IF (nyc == 1) THEN
  DO ic=1,ncx
    rhsc(ic,1) = 0.0
  END DO
END IF
IF (nyd == 1) THEN
  DO ic=1,ncx
    rhsc(ic,ncy) = 0.0
  END DO
END IF
RETURN
END SUBROUTINE res2

!     prolon2 modified from rgrd2u 11/20/97

SUBROUTINE prolon2(ncx,ncy,p,nx,ny,q,nxa,nxb,nyc,nyd,intpol)
IMPLICIT NONE
INTEGER, INTENT(IN OUT)                  :: ncx
INTEGER, INTENT(IN OUT)                      :: ncy
REAL (KIND=qPrec), INTENT(IN OUT)                     :: p(0:ncx+1,0:ncy+1)
INTEGER, INTENT(IN OUT)                      :: nx
INTEGER, INTENT(IN OUT)                  :: ny
REAL (KIND=qPrec), INTENT(IN OUT)                        :: q(0:nx+1,0:ny+1)
INTEGER, INTENT(IN OUT)                  :: nxa
INTEGER, INTENT(IN OUT)                  :: nxb
INTEGER, INTENT(IN OUT)                  :: nyc
INTEGER, INTENT(IN OUT)                  :: nyd
INTEGER, INTENT(IN OUT)                  :: intpol



INTEGER :: i,j,jc,ist,ifn,jst,jfn,joddst,joddfn

ist = 1
ifn = nx
jst = 1
jfn = ny
joddst = 1
joddfn = ny
IF (nxa == 1) THEN
  ist = 2
END IF
IF (nxb == 1) THEN
  ifn = nx-1
END IF
IF (nyc == 1) THEN
  jst = 2
  joddst = 3
END IF
IF (nyd == 1) THEN
  jfn = ny-1
  joddfn = ny-2
END IF
IF (intpol == 1 .OR. ncy < 4) THEN
  
!     linearly interpolate in y
  
  IF (ncy < ny) THEN
    
!     ncy grid is an every other point subset of ny grid
!     set odd j lines interpolating in x and then set even
!     j lines by averaging odd j lines
    
    DO j=joddst,joddfn,2
      jc = j/2+1
      CALL prolon1(ncx,p(0,jc),nx,q(0,j),nxa,nxb,intpol)
    END DO
    DO j=2,jfn,2
      DO i=ist,ifn
        q(i,j) = 0.5*(q(i,j-1)+q(i,j+1))
      END DO
    END DO
    
!     set periodic virtual boundaries if necessary
    
    IF (nyc == 0) THEN
      DO i=ist,ifn
        q(i,0) = q(i,ny-1)
        q(i,ny+1) = q(i,2)
      END DO
    END IF
    RETURN
  ELSE
    
!     ncy grid equals ny grid so interpolate in x only
    
    DO j=jst,jfn
      jc = j
      CALL prolon1(ncx,p(0,jc),nx,q(0,j),nxa,nxb,intpol)
    END DO
    
!     set periodic virtual boundaries if necessary
    
    IF (nyc == 0) THEN
      DO i=ist,ifn
        q(i,0) = q(i,ny-1)
        q(i,ny+1) = q(i,2)
      END DO
    END IF
    RETURN
  END IF
ELSE
  
!     cubically interpolate in y
  
  IF (ncy < ny) THEN
    
!     set every other point of ny grid by interpolating in x
    
    DO j=joddst,joddfn,2
      jc = j/2+1
      CALL prolon1(ncx,p(0,jc),nx,q(0,j),nxa,nxb,intpol)
    END DO
    
!     set deep interior of ny grid using values just
!     generated and symmetric cubic interpolation in y
    
    DO j=4,ny-3,2
      DO i=ist,ifn
        q(i,j)=(-q(i,j-3)+9.*(q(i,j-1)+q(i,j+1))-q(i,j+3))*.0625
      END DO
    END DO
    
!     interpolate from q at j=2 and j=ny-1
    
    IF (nyc /= 0) THEN
      
!     asymmetric formula near nonperiodic y boundaries
      
      DO i=ist,ifn
        q(i,2)=(5.*q(i,1)+15.*q(i,3)-5.*q(i,5)+q(i,7))*.0625
        q(i,ny-1)=(5.*q(i,ny)+15.*q(i,ny-2)-5.*q(i,ny-4)+ q(i,ny-6))*.0625
      END DO
    ELSE
      
!     periodicity in y alows symmetric formula near bndys
      
      DO i=ist,ifn
        q(i,2) = (-q(i,ny-2)+9.*(q(i,1)+q(i,3))-q(i,5))*.0625
        q(i,ny-1)=(-q(i,ny-4)+9.*(q(i,ny-2)+q(i,ny))-q(i,3))*.0625
        q(i,ny+1) = q(i,2)
        q(i,0) = q(i,ny-1)
      END DO
    END IF
    RETURN
  ELSE
    
!     ncy grid is equals ny grid so interpolate in x only
    
    DO j=jst,jfn
      jc = j
      CALL prolon1(ncx,p(0,jc),nx,q(0,j),nxa,nxb,intpol)
    END DO
    
!     set periodic virtual boundaries if necessary
    
    IF (nyc == 0) THEN
      DO i=ist,ifn
        q(i,0) = q(i,ny-1)
        q(i,ny+1) = q(i,2)
      END DO
    END IF
    RETURN
  END IF
END IF
END SUBROUTINE prolon2


!     11/20/97  modification of rgrd1u.f for mudpack

SUBROUTINE prolon1(ncx,p,nx,q,nxa,nxb,intpol)
IMPLICIT NONE
INTEGER, INTENT(IN OUT)                      :: ncx
REAL (KIND=qPrec), INTENT(IN OUT)                         :: p(0:ncx+1)
INTEGER, INTENT(IN OUT)                      :: nx
REAL (KIND=qPrec), INTENT(IN OUT)                        :: q(0:nx+1)
INTEGER, INTENT(IN OUT)                  :: nxa
INTEGER, INTENT(IN OUT)                  :: nxb
INTEGER, INTENT(IN OUT)                  :: intpol

INTEGER :: i,ic,ist,ifn,ioddst,ioddfn


ist = 1
ioddst = 1
ifn = nx
ioddfn = nx
IF (nxa == 1) THEN
  ist = 2
  ioddst = 3
END IF
IF (nxb == 1) THEN
  ifn = nx-1
  ioddfn = nx-2
END IF
IF (intpol == 1 .OR. ncx < 4) THEN
  
!     linear interpolation in x
  
  IF (ncx < nx) THEN
    
!     every other point of nx grid is ncx grid
    
    DO i=ioddst,ioddfn,2
      ic = (i+1)/2
      q(i) = p(ic)
    END DO
    DO i=2,ifn,2
      q(i) = 0.5*(q(i-1)+q(i+1))
    END DO
  ELSE
    
!     nx grid equals ncx grid
    
    DO i=ist,ifn
      q(i) = p(i)
    END DO
  END IF
  
!     set virtual end points if periodic
  
  IF (nxa == 0) THEN
    q(0) = q(nx-1)
    q(nx+1) = q(2)
  END IF
  RETURN
ELSE
  
!     cubic interpolation in x
  
  IF (ncx < nx) THEN
    DO i=ioddst,ioddfn,2
      ic = (i+1)/2
      q(i) = p(ic)
    END DO
    
!      set deep interior with symmetric formula
    
    DO i=4,nx-3,2
      q(i)=(-q(i-3)+9.*(q(i-1)+q(i+1))-q(i+3))*.0625
    END DO
    
!     interpolate from q at i=2 and i=nx-1
    
    IF (nxa /= 0) THEN
      
!     asymmetric formula near nonperiodic bndys
      
      q(2)=(5.*q(1)+15.*q(3)-5.*q(5)+q(7))*.0625
      q(nx-1)=(5.*q(nx)+15.*q(nx-2)-5.*q(nx-4)+q(nx-6))*.0625
    ELSE
      
!     periodicity in x alows symmetric formula near bndys
      
      q(2) = (-q(nx-2)+9.*(q(1)+q(3))-q(5))*.0625
      q(nx-1) = (-q(nx-4)+9.*(q(nx-2)+q(nx))-q(3))*.0625
      q(nx+1) = q(2)
      q(0) = q(nx-1)
    END IF
    RETURN
  ELSE
    
!     ncx grid equals nx grid
    
    DO i=ist,ifn
      q(i) = p(i)
    END DO
    IF (nxa == 0) THEN
      q(0) = q(nx-1)
      q(nx+1) = q(2)
    END IF
    RETURN
  END IF
END IF
END SUBROUTINE prolon1


SUBROUTINE cor2(nx,ny,phif,ncx,ncy,phic,nxa,nxb,nyc,nyd,intpol, phcor)

!     add coarse grid correction in phic to fine grid approximation
!     in phif using linear or cubic interpolation

IMPLICIT NONE
INTEGER, INTENT(IN OUT)                      :: nx
INTEGER, INTENT(IN OUT)                  :: ny
REAL (KIND=qPrec), INTENT(IN OUT)                        :: phif(0:nx+1,0:ny+1)
INTEGER, INTENT(IN OUT)                  :: ncx
INTEGER, INTENT(IN OUT)                  :: ncy
REAL (KIND=qPrec), INTENT(IN OUT)                     :: phic(0:ncx+1,0:ncy+1)
INTEGER, INTENT(IN OUT)                  :: nxa
INTEGER, INTENT(IN OUT)                  :: nxb
INTEGER, INTENT(IN OUT)                  :: nyc
INTEGER, INTENT(IN OUT)                  :: nyd
INTEGER, INTENT(IN OUT)                  :: intpol
REAL (KIND=qPrec), INTENT(IN OUT)                        :: phcor(0:nx+1,0:ny+1)

INTEGER :: i,j, ist,ifn,jst,jfn



DO j=0,ny+1
  DO i=0,nx+1
    phcor(i,j) = 0.0
  END DO
END DO

!     lift correction in phic to fine grid in phcor

CALL prolon2(ncx,ncy,phic,nx,ny,phcor,nxa,nxb,nyc,nyd,intpol)

!     add correction in phcor to phif on nonspecified boundaries

ist = 1
ifn = nx
jst = 1
jfn = ny
IF (nxa == 1) ist = 2
IF (nxb == 1) ifn = nx-1
IF (nyc == 1) jst = 2
IF (nyd == 1) jfn = ny-1
DO j=jst,jfn
  DO i=ist,ifn
    phif(i,j) = phif(i,j) + phcor(i,j)
  END DO
END DO

!     add periodic points if necessary

IF (nyc == 0) THEN
  DO i=ist,ifn
    phif(i,0) = phif(i,ny-1)
    phif(i,ny+1) = phif(i,2)
  END DO
END IF
IF (nxa == 0) THEN
  DO j=jst,jfn
    phif(0,j) = phif(nx-1,j)
    phif(nx+1,j) = phif(2,j)
  END DO
END IF
END SUBROUTINE cor2

SUBROUTINE pde2(nx,ny,u,i,j,ux3,ux4,uy3,uy4,nxa,nyc)
IMPLICIT NONE
INTEGER, INTENT(IN)                      :: nx
INTEGER, INTENT(IN)                      :: ny
REAL (KIND=qPrec), INTENT(IN OUT)                         :: u(nx,ny)
INTEGER, INTENT(IN OUT)                      :: i
INTEGER, INTENT(IN OUT)                      :: j
REAL (KIND=qPrec), INTENT(IN OUT)                        :: ux3
REAL (KIND=qPrec), INTENT(IN OUT)                        :: ux4
REAL (KIND=qPrec), INTENT(IN OUT)                        :: uy3
REAL (KIND=qPrec), INTENT(IN OUT)                        :: uy4
INTEGER, INTENT(IN OUT)                  :: nxa
INTEGER, INTENT(IN OUT)                  :: nyc


REAL (KIND=qPrec) :: dlx,dly,dlxx,dlyy,tdlx3,tdly3,dlx4,dly4
COMMON/pde2com/dlx,dly,dlxx,dlyy,tdlx3,tdly3,dlx4,dly4


!     use second order approximation in u to estimate (second order)
!     third and fourth partial derivatives in the x and y direction
!     non-symmetric difference formula (derived from the  routine
!     finpdf,findif) are used at and one point in from mixed boundaries.

IF (nxa /= 0) THEN
  
!     nonperiodic in x
  
  IF(i > 2 .AND. i < nx-1) THEN
    ux3 = (-u(i-2,j)+2.0*u(i-1,j)-2.0*u(i+1,j)+u(i+2,j))/tdlx3
    ux4 = (u(i-2,j)-4.0*u(i-1,j)+6.0*u(i,j)-4.0*u(i+1,j)+u(i+2,j)) /dlx4
  ELSE IF (i == 1) THEN
    ux3 = (-5.0*u(1,j)+18.0*u(2,j)-24.0*u(3,j)+14.0*u(4,j)- 3.0*u(5,j))/tdlx3
    ux4 = (3.0*u(1,j)-14.0*u(2,j)+26.0*u(3,j)-24.0*u(4,j)+  &
        11.0*u(5,j)-2.0*u(6,j))/dlx4
  ELSE IF (i == 2) THEN
    ux3 = (-3.0*u(1,j)+10.0*u(2,j)-12.0*u(3,j)+6.0*u(4,j)-u(5,j)) /tdlx3
    ux4 = (2.0*u(1,j)-9.0*u(2,j)+16.0*u(3,j)-14.0*u(4,j)+  &
        6.0*u(5,j)-u(6,j))/dlx4
  ELSE IF (i == nx-1) THEN
    ux3 = (u(nx-4,j)-6.0*u(nx-3,j)+12.0*u(nx-2,j)-10.0*u(nx-1,j)+  &
        3.0*u(nx,j))/tdlx3
    ux4 = (-u(nx-5,j)+6.0*u(nx-4,j)-14.0*u(nx-3,j)+16.0*u(nx-2,j)-  &
        9.0*u(nx-1,j)+2.0*u(nx,j))/dlx4
  ELSE IF (i == nx) THEN
    ux3 = (3.0*u(nx-4,j)-14.0*u(nx-3,j)+24.0*u(nx-2,j)-  &
        18.0*u(nx-1,j)+5.0*u(nx,j))/tdlx3
    ux4 = (-2.0*u(nx-5,j)+11.0*u(nx-4,j)-24.0*u(nx-3,j)+  &
        26.0*u(nx-2,j)-14.0*u(nx-1,j)+3.0*u(nx,j))/dlx4
  END IF
ELSE
  
!     periodic in x
  
  IF(i > 2 .AND. i < nx-1) THEN
    ux3 = (-u(i-2,j)+2.0*u(i-1,j)-2.0*u(i+1,j)+u(i+2,j))/tdlx3
    ux4 = (u(i-2,j)-4.0*u(i-1,j)+6.0*u(i,j)-4.0*u(i+1,j)+u(i+2,j)) /dlx4
  ELSE IF (i == 1) THEN
    ux3 = (-u(nx-2,j)+2.0*u(nx-1,j)-2.0*u(2,j)+u(3,j))/tdlx3
    ux4 = (u(nx-2,j)-4.0*u(nx-1,j)+6.0*u(1,j)-4.0*u(2,j)+u(3,j)) /dlx4
  ELSE IF (i == 2) THEN
    ux3 = (-u(nx-1,j)+2.0*u(1,j)-2.0*u(3,j)+u(4,j))/(tdlx3)
    ux4 = (u(nx-1,j)-4.0*u(1,j)+6.0*u(2,j)-4.0*u(3,j)+u(4,j))/dlx4
  ELSE IF (i == nx-1) THEN
    ux3 = (-u(nx-3,j)+2.0*u(nx-2,j)-2.0*u(1,j)+u(2,j))/tdlx3
    ux4 = (u(nx-3,j)-4.0*u(nx-2,j)+6.0*u(nx-1,j)-4.0*u(1,j)+ u(2,j))/dlx4
  ELSE IF (i == nx) THEN
    ux3 = (-u(nx-2,j)+2.0*u(nx-1,j)-2.0*u(2,j)+u(3,j))/tdlx3
    ux4 = (u(nx-2,j)-4.0*u(nx-1,j)+6.0*u(nx,j)-4.0*u(2,j)+u(3,j)) /dlx4
  END IF
END IF

!     y partial derivatives

IF (nyc /= 0) THEN
  
!     not periodic in y
  
  IF (j > 2 .AND. j < ny-1) THEN
    uy3 = (-u(i,j-2)+2.0*u(i,j-1)-2.0*u(i,j+1)+u(i,j+2))/tdly3
    uy4 = (u(i,j-2)-4.0*u(i,j-1)+6.0*u(i,j)-4.0*u(i,j+1)+u(i,j+2)) /dly4
  ELSE IF (j == 1) THEN
    uy3 = (-5.0*u(i,1)+18.0*u(i,2)-24.0*u(i,3)+14.0*u(i,4)- 3.0*u(i,5))/tdly3
    uy4 = (3.0*u(i,1)-14.0*u(i,2)+26.0*u(i,3)-24.0*u(i,4)+  &
        11.0*u(i,5)-2.0*u(i,6))/dly4
  ELSE IF (j == 2) THEN
    uy3 = (-3.0*u(i,1)+10.0*u(i,2)-12.0*u(i,3)+6.0*u(i,4)-u(i,5)) /tdly3
    uy4 = (2.0*u(i,1)-9.0*u(i,2)+16.0*u(i,3)-14.0*u(i,4)+  &
        6.0*u(i,5)-u(i,6))/dly4
  ELSE IF (j == ny-1) THEN
    uy3 = (u(i,ny-4)-6.0*u(i,ny-3)+12.0*u(i,ny-2)-10.0*u(i,ny-1)+  &
        3.0*u(i,ny))/tdly3
    uy4 = (-u(i,ny-5)+6.0*u(i,ny-4)-14.0*u(i,ny-3)+16.0*u(i,ny-2)-  &
        9.0*u(i,ny-1)+2.0*u(i,ny))/dly4
  ELSE IF (j == ny) THEN
    uy3 = (3.0*u(i,ny-4)-14.0*u(i,ny-3)+24.0*u(i,ny-2)-  &
        18.0*u(i,ny-1)+5.0*u(i,ny))/tdly3
    uy4 = (-2.0*u(i,ny-5)+11.0*u(i,ny-4)-24.0*u(i,ny-3)+  &
        26.0*u(i,ny-2)-14.0*u(i,ny-1)+3.0*u(i,ny))/dly4
  END IF
ELSE
  
!     periodic in y
  
  IF (j > 2 .AND. j < ny-1) THEN
    uy3 = (-u(i,j-2)+2.0*u(i,j-1)-2.0*u(i,j+1)+u(i,j+2))/tdly3
    uy4 = (u(i,j-2)-4.0*u(i,j-1)+6.0*u(i,j)-4.0*u(i,j+1)+u(i,j+2)) /dly4
  ELSE IF (j == 1) THEN
    uy3 = (-u(i,ny-2)+2.0*u(i,ny-1)-2.0*u(i,2)+u(i,3))/tdly3
    uy4 = (u(i,ny-2)-4.0*u(i,ny-1)+6.0*u(i,1)-4.0*u(i,2)+u(i,3)) /dly4
  ELSE IF (j == 2) THEN
    uy3 = (-u(i,ny-1)+2.0*u(i,1)-2.0*u(i,3)+u(i,4))/(tdly3)
    uy4 = (u(i,ny-1)-4.0*u(i,1)+6.0*u(i,2)-4.0*u(i,3)+u(i,4))/dly4
  ELSE IF (j == ny-1) THEN
    uy3 = (-u(i,ny-3)+2.0*u(i,ny-2)-2.0*u(i,1)+u(i,2))/tdly3
    uy4 = (u(i,ny-3)-4.0*u(i,ny-2)+6.0*u(i,ny-1)-4.0*u(i,1)+ u(i,2))/dly4
  ELSE IF (j == ny) THEN
    uy3 = (-u(i,ny-2)+2.0*u(i,ny-1)-2.0*u(i,2)+u(i,3))/tdly3
    uy4 = (u(i,ny-2)-4.0*u(i,ny-1)+6.0*u(i,ny)-4.0*u(i,2)+u(i,3)) /dly4
  END IF
END IF
RETURN
END SUBROUTINE pde2

SUBROUTINE swk3(nfx,nfy,nfz,phif,rhsf,phi,rhs)

!     set phif,rhsf input in arrays which include
!     virtual boundaries for phi (for all 2-d real codes)

IMPLICIT NONE
INTEGER, INTENT(IN OUT)                      :: nfx
INTEGER, INTENT(IN OUT)                  :: nfy
INTEGER, INTENT(IN OUT)                  :: nfz
REAL (KIND=qPrec), INTENT(IN OUT)                         :: phif(nfx,nfy,nfz)
REAL (KIND=qPrec), INTENT(IN OUT)                         :: rhsf(nfx,nfy,nfz)
REAL (KIND=qPrec), INTENT(IN OUT)                        :: phi(0:nfx+1,0:nfy+1,0:nfz+1)
REAL (KIND=qPrec), INTENT(IN OUT)                        :: rhs(nfx,nfy,nfz)

INTEGER :: i,j,k



DO k=1,nfz
  DO j=1,nfy
    DO i=1,nfx
      phi(i,j,k) = phif(i,j,k)
      rhs(i,j,k) = rhsf(i,j,k)
    END DO
  END DO
END DO

!     set virtual boundaries in phi to zero

DO k=0,nfz+1
  DO j=0,nfy+1
    phi(0,j,k) = 0.0
    phi(nfx+1,j,k) = 0.0
  END DO
END DO
DO k=0,nfz+1
  DO i=0,nfx+1
    phi(i,0,k) = 0.0
    phi(i,nfy+1,k) = 0.0
  END DO
END DO
DO j=0,nfy+1
  DO i=0,nfx+1
    phi(i,j,0) = 0.0
    phi(i,j,nfz+1) = 0.0
  END DO
END DO
RETURN
END SUBROUTINE swk3

SUBROUTINE trsfc3(nx,ny,nz,phi,rhs,ncx,ncy,ncz,phic,rhsc)

!     transfer fine grid to coarse grid

IMPLICIT NONE
INTEGER, INTENT(IN OUT)                      :: nx
INTEGER, INTENT(IN OUT)                      :: ny
INTEGER, INTENT(IN OUT)                  :: nz
REAL (KIND=qPrec), INTENT(IN OUT)                         :: phi(0:nx+1,0:ny+1,0:nz+1)
REAL (KIND=qPrec), INTENT(IN OUT)                         :: rhs(nx,ny,nz)
INTEGER, INTENT(IN OUT)                      :: ncx
INTEGER, INTENT(IN OUT)                  :: ncy
INTEGER, INTENT(IN OUT)                  :: ncz
REAL (KIND=qPrec), INTENT(IN OUT)                        :: phic(0:ncx+1,0:ncy+1,0:ncz+1)
REAL (KIND=qPrec), INTENT(IN OUT)                        :: rhsc(ncx,ncy,ncz)

INTEGER :: i,j,k,ic,jc,kc,ix,jy,kz



!     set virtual boundaries in phic to zero

DO kc=0,ncz+1
  DO jc=0,ncy+1
    phic(0,jc,kc) = 0.0
    phic(ncx+1,jc,kc) = 0.0
  END DO
END DO
DO kc=0,ncz+1
  DO ic=0,ncx+1
    phic(ic,0,kc) = 0.0
    phic(ic,ncy+1,kc) = 0.0
  END DO
END DO
DO jc=0,ncy+1
  DO ic=0,ncx+1
    phic(ic,jc,0) = 0.0
    phic(ic,jc,ncz+1) = 0.0
  END DO
END DO
IF (ncx < nx .AND. ncy < ny .AND. ncz < nz) THEN
  
!     coarsening in x,y,z (usually the case?)
  
  DO kc=1,ncz
    k = kc+kc-1
    DO jc=1,ncy
      j = jc+jc-1
      DO ic=1,ncx
        i = ic+ic-1
        phic(ic,jc,kc) = phi(i,j,k)
        rhsc(ic,jc,kc) = rhs(i,j,k)
      END DO
    END DO
  END DO
ELSE
  
!     no coarsening in at least one dimension
  
  ix = 1
  IF (ncx == nx) ix = 0
  jy = 1
  IF (ncy == ny) jy = 0
  kz = 1
  IF (ncz == nz) kz = 0
  
  DO kc=1,ncz
    k = kc+kz*(kc-1)
    DO jc=1,ncy
      j = jc+jy*(jc-1)
      DO ic=1,ncx
        i = ic+ix*(ic-1)
        phic(ic,jc,kc) = phi(i,j,k)
        rhsc(ic,jc,kc) = rhs(i,j,k)
      END DO
    END DO
  END DO
END IF
RETURN
END SUBROUTINE trsfc3

SUBROUTINE res3(nx,ny,nz,resf,ncx,ncy,ncz,rhsc, nxa,nxb,nyc,nyd,nze,nzf)
IMPLICIT NONE
INTEGER, INTENT(IN)                      :: nx
INTEGER, INTENT(IN)                      :: ny
INTEGER, INTENT(IN)                      :: nz
REAL (KIND=qPrec), INTENT(IN OUT)                         :: resf(nx,ny,nz)
INTEGER, INTENT(IN)                      :: ncx
INTEGER, INTENT(IN)                      :: ncy
INTEGER, INTENT(IN)                      :: ncz
REAL (KIND=qPrec), INTENT(IN OUT)                        :: rhsc(ncx,ncy,ncz)
INTEGER, INTENT(IN OUT)                      :: nxa
INTEGER, INTENT(IN OUT)                      :: nxb
INTEGER, INTENT(IN OUT)                      :: nyc
INTEGER, INTENT(IN OUT)                      :: nyd
INTEGER, INTENT(IN OUT)                      :: nze
INTEGER, INTENT(IN OUT)                      :: nzf


INTEGER :: ix,jy,kz,i,j,k,ic,jc,kc,im1,ip1,jm1,jp1,km1,kp1
REAL (KIND=qPrec) :: rm,rk,rp

!     restrict fine grid residual in resf to coarse grid in rhsc
!     using full weighting



!     set x,y,z coarsening integer subscript scales

ix = 1
IF (ncx == nx) ix = 0
jy = 1
IF (ncy == ny) jy = 0
kz = 1
IF (ncz == nz) kz = 0

!     restrict on interior

IF (ncz < nz .AND. ncy < ny .AND. ncx < nx) THEN
  
!     coarsening in x,y,z
  
!$OMP PARALLEL DO PRIVATE(i,j,k,ic,jc,kc,rm,rk,rp)
!$OMP+SHARED(resf,rhsc,ncx,ncy,ncz)
  DO kc=2,ncz-1
    k = kc+kc-1
    DO jc=2,ncy-1
      j = jc+jc-1
      DO ic=2,ncx-1
        i = ic+ic-1
        
!     weight on k-1,k,k+1 z planes in rm,rk,rp
        
        rm=(resf(i-1,j-1,k-1)+resf(i+1,j-1,k-1)+resf(i-1,j+1,k-1)+  &
            resf(i+1,j+1,k-1)+2.*(resf(i-1,j,k-1)+resf(i+1,j,k-1)+  &
            resf(i,j-1,k-1)+resf(i,j+1,k-1))+4.*resf(i,j,k-1))*.0625
        rk=(resf(i-1,j-1,k)+resf(i+1,j-1,k)+resf(i-1,j+1,k)+  &
            resf(i+1,j+1,k)+2.*(resf(i-1,j,k)+resf(i+1,j,k)+  &
            resf(i,j-1,k)+resf(i,j+1,k))+4.*resf(i,j,k))*.0625
        rp=(resf(i-1,j-1,k+1)+resf(i+1,j-1,k+1)+resf(i-1,j+1,k+1)+  &
            resf(i+1,j+1,k+1)+2.*(resf(i-1,j,k+1)+resf(i+1,j,k+1)+  &
            resf(i,j-1,k+1)+resf(i,j+1,k+1))+4.*resf(i,j,k+1))*.0625
        
!     weight in z direction for final result
        
        rhsc(ic,jc,kc) = 0.25*(rm+2.*rk+rp)
      END DO
    END DO
  END DO
ELSE
  
!     allow for noncoarsening in any of x,y,z
  
!$OMP PARALLEL DO PRIVATE(i,j,k,ic,jc,kc,rm,rk,rp)
!$OMP+SHARED(ix,jy,kz,resf,rhsc,ncx,ncy,ncz)
  DO kc=2,ncz-1
    k = kc+kz*(kc-1)
    DO jc=2,ncy-1
      j = jc+jy*(jc-1)
      DO ic=2,ncx-1
        i = ic+ix*(ic-1)
        
!     weight on k-1,k,k+1 z planes in rm,rk,rp
        
        rm=(resf(i-1,j-1,k-1)+resf(i+1,j-1,k-1)+resf(i-1,j+1,k-1)+  &
            resf(i+1,j+1,k-1)+2.*(resf(i-1,j,k-1)+resf(i+1,j,k-1)+  &
            resf(i,j-1,k-1)+resf(i,j+1,k-1))+4.*resf(i,j,k-1))*.0625
        rk=(resf(i-1,j-1,k)+resf(i+1,j-1,k)+resf(i-1,j+1,k)+  &
            resf(i+1,j+1,k)+2.*(resf(i-1,j,k)+resf(i+1,j,k)+  &
            resf(i,j-1,k)+resf(i,j+1,k))+4.*resf(i,j,k))*.0625
        rp=(resf(i-1,j-1,k+1)+resf(i+1,j-1,k+1)+resf(i-1,j+1,k+1)+  &
            resf(i+1,j+1,k+1)+2.*(resf(i-1,j,k+1)+resf(i+1,j,k+1)+  &
            resf(i,j-1,k+1)+resf(i,j+1,k+1))+4.*resf(i,j,k+1))*.0625
        
!     weight in z direction for final result
        
        rhsc(ic,jc,kc) = 0.25*(rm+2.*rk+rp)
      END DO
    END DO
  END DO
END IF

!     set residual on boundaries

DO ic=1,ncx,ncx-1
  
!     x=xa and x=xb
  
  i = ic+ix*(ic-1)
  im1 = MAX0(i-1,2)
  ip1 = MIN0(i+1,nx-1)
  IF (i == 1 .AND. nxa == 0) im1 = nx-1
  IF (i == nx .AND. nxb == 0) ip1 = 2
  
!    (y,z) interior
  
!$OMP PARALLEL DO PRIVATE(j,k,jc,kc,rm,rk,rp)
!$OMP+SHARED(kz,jy,ic,im1,i,ip1,resf,rhsc,ncy,ncz)
  DO kc=2,ncz-1
    k = kc+kz*(kc-1)
    DO jc=2,ncy-1
      j = jc+jy*(jc-1)
      rm=(resf(im1,j-1,k-1)+resf(ip1,j-1,k-1)+resf(im1,j+1,k-1)+  &
          resf(ip1,j+1,k-1)+2.*(resf(im1,j,k-1)+resf(ip1,j,k-1)+  &
          resf(i,j-1,k-1)+resf(i,j+1,k-1))+4.*resf(i,j,k-1))*.0625
      rk=(resf(im1,j-1,k)+resf(ip1,j-1,k)+resf(im1,j+1,k)+  &
          resf(ip1,j+1,k)+2.*(resf(im1,j,k)+resf(ip1,j,k)+  &
          resf(i,j-1,k)+resf(i,j+1,k))+4.*resf(i,j,k))*.0625
      rp=(resf(im1,j-1,k+1)+resf(ip1,j-1,k+1)+resf(im1,j+1,k+1)+  &
          resf(ip1,j+1,k+1)+2.*(resf(im1,j,k+1)+resf(ip1,j,k+1)+  &
          resf(i,j-1,k+1)+resf(i,j+1,k+1))+4.*resf(i,j,k+1))*.0625
      rhsc(ic,jc,kc) = 0.25*(rm+2.*rk+rp)
    END DO
  END DO
  
!     x=xa,xb and y=yc,yd interior edges
  
  DO jc=1,ncy,ncy-1
    j = jc+jy*(jc-1)
    jm1 = MAX0(j-1,2)
    jp1 = MIN0(j+1,ny-1)
    IF (j == 1 .AND. nyc == 0) jm1 = ny-1
    IF (j == ny .AND. nyc == 0) jp1 = 2
    DO kc=2,ncz-1
      k = kc+kz*(kc-1)
      rm=(resf(im1,jm1,k-1)+resf(ip1,jm1,k-1)+resf(im1,jp1,k-1)+  &
          resf(ip1,jp1,k-1)+2.*(resf(im1,j,k-1)+resf(ip1,j,k-1)+  &
          resf(i,jm1,k-1)+resf(i,jp1,k-1))+4.*resf(i,j,k-1))*.0625
      rk=(resf(im1,jm1,k)+resf(ip1,jm1,k)+resf(im1,jp1,k)+  &
          resf(ip1,jp1,k)+2.*(resf(im1,j,k)+resf(ip1,j,k)+  &
          resf(i,jm1,k)+resf(i,jp1,k))+4.*resf(i,j,k))*.0625
      rp=(resf(im1,jm1,k+1)+resf(ip1,jm1,k+1)+resf(im1,jp1,k+1)+  &
          resf(ip1,jp1,k+1)+2.*(resf(im1,j,k+1)+resf(ip1,j,k+1)+  &
          resf(i,jm1,k+1)+resf(i,jp1,k+1))+4.*resf(i,j,k+1))*.0625
      rhsc(ic,jc,kc) = 0.25*(rm+2.*rk+rp)
    END DO
!     x=xa,xb; y=yc,yd; z=ze,zf cornors
    DO kc=1,ncz,ncz-1
      k = kc+kz*(kc-1)
      km1 = MAX0(k-1,2)
      kp1 = MIN0(k+1,nz-1)
      IF (k == 1 .AND. nze == 0) km1 = nz-1
      IF (k == nz .AND. nzf == 0) kp1 = 2
      rm=(resf(im1,jm1,km1)+resf(ip1,jm1,km1)+resf(im1,jp1,km1)+  &
          resf(ip1,jp1,km1)+2.*(resf(im1,j,km1)+resf(ip1,j,km1)+  &
          resf(i,jm1,km1)+resf(i,jp1,km1))+4.*resf(i,j,km1))*.0625
      rk=(resf(im1,jm1,k)+resf(ip1,jm1,k)+resf(im1,jp1,k)+  &
          resf(ip1,jp1,k)+2.*(resf(im1,j,k)+resf(ip1,j,k)+  &
          resf(i,jm1,k)+resf(i,jp1,k))+4.*resf(i,j,k))*.0625
      rp=(resf(im1,jm1,kp1)+resf(ip1,jm1,kp1)+resf(im1,jp1,kp1)+  &
          resf(ip1,jp1,kp1)+2.*(resf(im1,j,kp1)+resf(ip1,j,kp1)+  &
          resf(i,jm1,kp1)+resf(i,jp1,kp1))+4.*resf(i,j,kp1))*.0625
      rhsc(ic,jc,kc) = 0.25*(rm+2.*rk+rp)
    END DO
  END DO
  
!      x=xa,xb and z=ze,zf edges
  
  DO kc=1,ncz,ncz-1
    k = kc+kz*(kc-1)
    km1 = MAX0(k-1,2)
    kp1 = MIN0(k+1,nz-1)
    IF (k == 1 .AND. nze == 0) km1 = nz-1
    IF (k == nz .AND. nzf == 0) kp1 = 2
    DO jc=2,ncy-1
      j = jc+jy*(jc-1)
      rm=(resf(im1,j-1,km1)+resf(ip1,j-1,km1)+resf(im1,j+1,km1)+  &
          resf(ip1,j+1,km1)+2.*(resf(im1,j,km1)+resf(ip1,j,km1)+  &
          resf(i,j-1,km1)+resf(i,j+1,km1))+4.*resf(i,j,km1))*.0625
      rk=(resf(im1,j-1,k)+resf(ip1,j-1,k)+resf(im1,j+1,k)+  &
          resf(ip1,j+1,k)+2.*(resf(im1,j,k)+resf(ip1,j,k)+  &
          resf(i,j-1,k)+resf(i,j+1,k))+4.*resf(i,j,k))*.0625
      rp=(resf(im1,j-1,kp1)+resf(ip1,j-1,kp1)+resf(im1,j+1,kp1)+  &
          resf(ip1,j+1,kp1)+2.*(resf(im1,j,kp1)+resf(ip1,j,kp1)+  &
          resf(i,j-1,kp1)+resf(i,j+1,kp1))+4.*resf(i,j,kp1))*.0625
      rhsc(ic,jc,kc) = 0.25*(rm+2.*rk+rp)
    END DO
  END DO
END DO

!     y boundaries y=yc and y=yd

DO jc=1,ncy,ncy-1
  j = jc+jy*(jc-1)
  jm1 = MAX0(j-1,2)
  jp1 = MIN0(j+1,ny-1)
  IF (j == 1 .AND. nyc == 0) jm1 = ny-1
  IF (j == ny .AND. nyd == 0) jp1 = 2
  
!     (x,z) interior
  
!$OMP PARALLEL DO PRIVATE(i,k,ic,kc,rm,rk,rp)
!$OMP+SHARED(ix,kz,jc,jm1,j,jp1,resf,rhsc,ncx,ncz)
  DO kc=2,ncz-1
    k = kc+kz*(kc-1)
    DO ic=2,ncx-1
      i = ic+ix*(ic-1)
      rm=(resf(i-1,jm1,k-1)+resf(i+1,jm1,k-1)+resf(i-1,jp1,k-1)+  &
          resf(i+1,jp1,k-1)+2.*(resf(i-1,j,k-1)+resf(i+1,j,k-1)+  &
          resf(i,jm1,k-1)+resf(i,jp1,k-1))+4.*resf(i,j,k-1))*.0625
      rk=(resf(i-1,jm1,k)+resf(i+1,jm1,k)+resf(i-1,jp1,k)+  &
          resf(i+1,jp1,k)+2.*(resf(i-1,j,k)+resf(i+1,j,k)+  &
          resf(i,jm1,k)+resf(i,jp1,k))+4.*resf(i,j,k))*.0625
      rp=(resf(i-1,jm1,k+1)+resf(i+1,jm1,k+1)+resf(i-1,jp1,k+1)+  &
          resf(i+1,jp1,k+1)+2.*(resf(i-1,j,k+1)+resf(i+1,j,k+1)+  &
          resf(i,jm1,k+1)+resf(i,jp1,k+1))+4.*resf(i,j,k+1))*.0625
      rhsc(ic,jc,kc) = 0.25*(rm+2.*rk+rp)
    END DO
  END DO
  
!     y=yc,yd and z=ze,zf edges
  
  DO kc=1,ncz,ncz-1
    k = kc+kz*(kc-1)
    km1 = MAX0(k-1,2)
    kp1 = MIN0(k+1,nz-1)
    IF (k == 1 .AND. nze == 0) km1 = nz-1
    IF (k == nz .AND. nzf == 0) kp1 = 2
    
!     interior in x
    
    DO ic=2,ncx-1
      i = ic+ix*(ic-1)
      rm=(resf(i-1,jm1,km1)+resf(i+1,jm1,km1)+resf(i-1,jp1,km1)+  &
          resf(i+1,jp1,km1)+2.*(resf(i-1,j,km1)+resf(i+1,j,km1)+  &
          resf(i,jm1,km1)+resf(i,jp1,km1))+4.*resf(i,j,km1))*.0625
      rk=(resf(i-1,jm1,k)+resf(i+1,jm1,k)+resf(i-1,jp1,k)+  &
          resf(i+1,jp1,k)+2.*(resf(i-1,j,k)+resf(i+1,j,k)+  &
          resf(i,jm1,k)+resf(i,jp1,k))+4.*resf(i,j,k))*.0625
      rp=(resf(i-1,jm1,kp1)+resf(i+1,jm1,kp1)+resf(i-1,jp1,kp1)+  &
          resf(i+1,jp1,kp1)+2.*(resf(i-1,j,kp1)+resf(i+1,j,kp1)+  &
          resf(i,jm1,kp1)+resf(i,jp1,kp1))+4.*resf(i,j,kp1))*.0625
      rhsc(ic,jc,kc) = 0.25*(rm+2.*rk+rp)
    END DO
  END DO
END DO

!     z=ze,zf boundaries

DO kc=1,ncz,ncz-1
  k = kc+kz*(kc-1)
  km1 = MAX0(k-1,2)
  kp1 = MIN0(k+1,nz-1)
  IF (k == 1 .AND. nze == 0) km1 = nz-1
  IF (k == nz .AND. nzf == 0) kp1 = 2
  
!     (x,y) interior
  
!$OMP PARALLEL DO PRIVATE(i,j,ic,jc,rm,rk,rp)
!$OMP+SHARED(ix,jy,kc,km1,k,kp1,resf,rhsc,ncx,ncz)
  DO jc=2,ncy-1
    j = jc+jy*(jc-1)
    DO ic=2,ncx-1
      i = ic+ix*(ic-1)
      rm=(resf(i-1,j-1,km1)+resf(i+1,j-1,km1)+resf(i-1,j+1,km1)+  &
          resf(i+1,j+1,km1)+2.*(resf(i-1,j,km1)+resf(i+1,j,km1)+  &
          resf(i,j-1,km1)+resf(i,j+1,km1))+4.*resf(i,j,km1))*.0625
      rk=(resf(i-1,j-1,k)+resf(i+1,j-1,k)+resf(i-1,j+1,k)+  &
          resf(i+1,j+1,k)+2.*(resf(i-1,j,k)+resf(i+1,j,k)+  &
          resf(i,j-1,k)+resf(i,j+1,k))+4.*resf(i,j,k))*.0625
      rp=(resf(i-1,j-1,kp1)+resf(i+1,j-1,kp1)+resf(i-1,j+1,kp1)+  &
          resf(i+1,j+1,kp1)+2.*(resf(i-1,j,kp1)+resf(i+1,j,kp1)+  &
          resf(i,j-1,kp1)+resf(i,j+1,kp1))+4.*resf(i,j,kp1))*.0625
      rhsc(ic,jc,kc) = 0.25*(rm+2.*rk+rp)
    END DO
  END DO
END DO

!     set coarse grid residual to zero at specified boundaries

IF (nxa == 1) THEN
  ic = 1
  DO kc=1,ncz
    DO jc=1,ncy
      rhsc(ic,jc,kc) = 0.0
    END DO
  END DO
END IF
IF (nxb == 1) THEN
  ic = ncx
  DO kc=1,ncz
    DO jc=1,ncy
      rhsc(ic,jc,kc) = 0.0
    END DO
  END DO
END IF
IF (nyc == 1) THEN
  jc = 1
  DO kc=1,ncz
    DO ic=1,ncx
      rhsc(ic,jc,kc) = 0.0
    END DO
  END DO
END IF
IF (nyd == 1) THEN
  jc = ncy
  DO kc=1,ncz
    DO ic=1,ncx
      rhsc(ic,jc,kc) = 0.0
    END DO
  END DO
END IF
IF (nze == 1) THEN
  kc = 1
  DO jc=1,ncy
    DO ic=1,ncx
      rhsc(ic,jc,kc) = 0.0
    END DO
  END DO
END IF
IF (nzf == 1) THEN
  kc = ncz
  DO jc=1,ncy
    DO ic=1,ncx
      rhsc(ic,jc,kc) = 0.0
    END DO
  END DO
END IF
RETURN
END SUBROUTINE res3


!     prolon3 modified from prolon2 11/25/97

SUBROUTINE prolon3(ncx,ncy,ncz,p,nx,ny,nz,q,nxa,nxb,nyc,nyd, nze,nzf,intpol)
IMPLICIT NONE
INTEGER, INTENT(IN OUT)                  :: ncx
INTEGER, INTENT(IN OUT)                  :: ncy
INTEGER, INTENT(IN OUT)                      :: ncz
REAL (KIND=qPrec), INTENT(IN OUT)                     :: p(0:ncx+1,0:ncy+1,0:ncz+1)
INTEGER, INTENT(IN OUT)                      :: nx
INTEGER, INTENT(IN OUT)                      :: ny
INTEGER, INTENT(IN OUT)                  :: nz
REAL (KIND=qPrec), INTENT(IN OUT)                        :: q(0:nx+1,0:ny+1,0:nz+1)
INTEGER, INTENT(IN OUT)                  :: nxa
INTEGER, INTENT(IN OUT)                  :: nxb
INTEGER, INTENT(IN OUT)                  :: nyc
INTEGER, INTENT(IN OUT)                  :: nyd
INTEGER, INTENT(IN OUT)                  :: nze
INTEGER, INTENT(IN OUT)                  :: nzf
INTEGER, INTENT(IN OUT)                  :: intpol



INTEGER :: i,j,k,kc,ist,ifn,jst,jfn,kst,kfn,koddst,koddfn

ist = 1
ifn = nx
jst = 1
jfn = ny
kst = 1
kfn = nz
koddst = 1
koddfn = nz
IF (nxa == 1) THEN
  ist = 2
END IF
IF (nxb == 1) THEN
  ifn = nx-1
END IF
IF (nyc == 1) THEN
  jst = 2
END IF
IF (nyd == 1) THEN
  jfn = ny-1
END IF
IF (nze == 1) THEN
  kst = 2
  koddst = 3
END IF
IF (nzf == 1) THEN
  kfn = nz-1
  koddfn = nz-2
END IF
IF (intpol == 1 .OR. ncz < 4) THEN
  
!     linearly interpolate in z
  
  IF (ncz < nz) THEN
    
!     ncz grid is an every other point subset of nz grid
!     set odd k planes interpolating in x&y and then set even
!     k planes by averaging odd k planes
    
    DO k=koddst,koddfn,2
      kc = k/2+1
      CALL prolon2(ncx,ncy,p(0,0,kc),nx,ny,q(0,0,k),nxa,nxb,nyc, nyd,intpol)
    END DO
    DO k=2,kfn,2
      DO j=jst,jfn
        DO i=ist,ifn
          q(i,j,k) = 0.5*(q(i,j,k-1)+q(i,j,k+1))
        END DO
      END DO
    END DO
    
!     set periodic virtual boundaries if necessary
    
    IF (nze == 0) THEN
      DO j=jst,jfn
        DO i=ist,ifn
          q(i,j,0) = q(i,j,nz-1)
          q(i,j,nz+1) = q(i,j,2)
        END DO
      END DO
    END IF
    RETURN
  ELSE
    
!     ncz grid is equals nz grid so interpolate in x&y only
    
    DO k=kst,kfn
      kc = k
      CALL prolon2(ncx,ncy,p(0,0,kc),nx,ny,q(0,0,k),nxa,nxb,nyc, nyd,intpol)
    END DO
    
!     set periodic virtual boundaries if necessary
    
    IF (nze == 0) THEN
      DO j=jst,jfn
        DO i=ist,ifn
          q(i,j,0) = q(i,j,nz-1)
          q(i,j,nz+1) = q(i,j,2)
        END DO
      END DO
    END IF
    RETURN
  END IF
ELSE
  
!     cubically interpolate in z
  
  IF (ncz < nz) THEN
    
!     set every other point of nz grid by interpolating in x&y
    
    DO k=koddst,koddfn,2
      kc = k/2+1
      CALL prolon2(ncx,ncy,p(0,0,kc),nx,ny,q(0,0,k),nxa,nxb,nyc, nyd,intpol)
    END DO
    
!     set deep interior of nz grid using values just
!     generated and symmetric cubic interpolation in z
    
    DO k=4,nz-3,2
      DO j=jst,jfn
        DO i=ist,ifn
          q(i,j,k)=(-q(i,j,k-3)+9.*(q(i,j,k-1)+q(i,j,k+1))-q(i,j,k+3)) *.0625
        END DO
      END DO
    END DO
    
!     interpolate from q at k=2 and k=nz-1
    
    IF (nze /= 0) THEN
      
!     asymmetric formula near nonperiodic z boundaries
      
      DO j=jst,jfn
        DO i=ist,ifn
          q(i,j,2)=(5.*q(i,j,1)+15.*q(i,j,3)-5.*q(i,j,5)+q(i,j,7)) *.0625
          q(i,j,nz-1)=(5.*q(i,j,nz)+15.*q(i,j,nz-2)-5.*q(i,j,nz-4)+  &
              q(i,j,nz-6))*.0625
        END DO
      END DO
    ELSE
      
!     periodicity in y alows symmetric formula near bndys
      
      DO j=jst,jfn
        DO i=ist,ifn
          q(i,j,2) = (-q(i,j,nz-2)+9.*(q(i,j,1)+q(i,j,3))-q(i,j,5)) *.0625
          q(i,j,nz-1)=(-q(i,j,nz-4)+9.*(q(i,j,nz-2)+q(i,j,nz))-  &
              q(i,j,3))*.0625
          q(i,j,nz+1) = q(i,j,2)
          q(i,j,0) = q(i,j,nz-1)
        END DO
      END DO
    END IF
    RETURN
  ELSE
    
!     ncz grid is equals nx grid so interpolate in x&y only
    
    DO k=kst,kfn
      kc = k
      CALL prolon2(ncx,ncy,p(0,0,kc),nx,ny,q(0,0,k),nxa,nxb,nyc, nyd,intpol)
    END DO
    
!     set periodic virtual boundaries if necessary
    
    IF (nze == 0) THEN
      DO j=jst,jfn
        DO i=ist,ifn
          q(i,j,0) = q(i,j,nz-1)
          q(i,j,nz+1) = q(i,j,2)
        END DO
      END DO
    END IF
    RETURN
  END IF
END IF
END SUBROUTINE prolon3

SUBROUTINE cor3(nx,ny,nz,phif,ncx,ncy,ncz,phic,nxa,nxb,nyc,nyd,  &
    nze,nzf,intpol,phcor)
IMPLICIT NONE
INTEGER, INTENT(IN OUT)                      :: nx
INTEGER, INTENT(IN OUT)                  :: ny
INTEGER, INTENT(IN OUT)                  :: nz
REAL (KIND=qPrec), INTENT(IN OUT)                        :: phif(0:nx+1,0:ny+1,0:nz+1)
INTEGER, INTENT(IN OUT)                  :: ncx
INTEGER, INTENT(IN OUT)                  :: ncy
INTEGER, INTENT(IN OUT)                  :: ncz
REAL (KIND=qPrec), INTENT(IN OUT)                     :: phic(0:ncx+1,0:ncy+1,0:ncz+1)
INTEGER, INTENT(IN OUT)                  :: nxa
INTEGER, INTENT(IN OUT)                  :: nxb
INTEGER, INTENT(IN OUT)                  :: nyc
INTEGER, INTENT(IN OUT)                  :: nyd
INTEGER, INTENT(IN OUT)                  :: nze
INTEGER, INTENT(IN OUT)                  :: nzf
INTEGER, INTENT(IN OUT)                  :: intpol
REAL (KIND=qPrec), INTENT(IN OUT)                        :: phcor(0:nx+1,0:ny+1,0:nz+1)


INTEGER :: i,j,k,ist,ifn,jst,jfn,kst,kfn

!     add coarse grid correction in phic to fine grid approximation
!     in phif using linear or cubic interpolation




DO k=0,nz+1
  DO j=0,ny+1
    DO i=0,nx+1
      phcor(i,j,k) = 0.0
    END DO
  END DO
END DO

!     lift correction in phic to fine grid in phcor

CALL prolon3(ncx,ncy,ncz,phic,nx,ny,nz,phcor,nxa,nxb,nyc,nyd, nze,nzf,intpol)

!     add correction in phcor to phif on nonspecified boundaries

ist = 1
ifn = nx
jst = 1
jfn = ny
kst = 1
kfn = nz
IF (nxa == 1) ist = 2
IF (nxb == 1) ifn = nx-1
IF (nyc == 1) jst = 2
IF (nyd == 1) jfn = ny-1
IF (nze == 1) kst = 2
IF (nzf == 1) kfn = nz-1
DO k=kst,kfn
  DO j=jst,jfn
    DO i=ist,ifn
      phif(i,j,k) = phif(i,j,k) + phcor(i,j,k)
    END DO
  END DO
END DO

!     add periodic points if necessary

IF (nze == 0) THEN
  DO j=jst,jfn
    DO i=ist,ifn
      phif(i,j,0) = phif(i,j,nz-1)
      phif(i,j,nz+1) = phif(i,j,2)
    END DO
  END DO
END IF
IF (nyc == 0) THEN
  DO k=kst,kfn
    DO i=ist,ifn
      phif(i,0,k) = phif(i,ny-1,k)
      phif(i,ny+1,k) = phif(i,2,k)
    END DO
  END DO
END IF
IF (nxa == 0) THEN
  DO k=kst,kfn
    DO j=jst,jfn
      phif(0,j,k) = phif(nx-1,j,k)
      phif(nx+1,j,k) = phif(2,j,k)
    END DO
  END DO
END IF
END SUBROUTINE cor3

SUBROUTINE per3vb(nx,ny,nz,phi,nxa,nyc,nze)

!     set virtual periodic boundaries from interior values
!     in three dimensions (for all 3-d solvers)

IMPLICIT NONE
INTEGER, INTENT(IN)                      :: nx
INTEGER, INTENT(IN)                  :: ny
INTEGER, INTENT(IN)                  :: nz
REAL (KIND=qPrec), INTENT(IN OUT)                        :: phi(0:nx+1,0:ny+1,0:nz+1)
INTEGER, INTENT(IN)                  :: nxa
INTEGER, INTENT(IN)                  :: nyc
INTEGER, INTENT(IN)                  :: nze

INTEGER :: j,k,i


IF (nxa == 0) THEN
  DO k=1,nz
    DO j=1,ny
      phi(0,j,k) = phi(nx-1,j,k)
      phi(nx,j,k) = phi(1,j,k)
      phi(nx+1,j,k) = phi(2,j,k)
    END DO
  END DO
END IF
IF (nyc == 0) THEN
  DO k=1,nz
    DO i=1,nx
      phi(i,0,k) = phi(i,ny-1,k)
      phi(i,ny,k) = phi(i,1,k)
      phi(i,ny+1,k) = phi(i,2,k)
    END DO
  END DO
END IF
IF (nze == 0) THEN
  DO j=1,ny
    DO i=1,nx
      phi(i,j,0) = phi(i,j,nz-1)
      phi(i,j,nz) = phi(i,j,1)
      phi(i,j,nz+1) = phi(i,j,2)
    END DO
  END DO
END IF
RETURN
END SUBROUTINE per3vb

SUBROUTINE pde2cr(nx,ny,u,i,j,ux3y,uxy3,ux2y2)

!     compute mixed partial derivative approximations

IMPLICIT NONE
INTEGER, INTENT(IN)                      :: nx
INTEGER, INTENT(IN)                      :: ny
REAL (KIND=qPrec), INTENT(IN OUT)                         :: u(nx,ny)
INTEGER, INTENT(IN OUT)                  :: i
INTEGER, INTENT(IN OUT)                      :: j
REAL (KIND=qPrec), INTENT(IN OUT)                        :: ux3y
REAL (KIND=qPrec), INTENT(IN OUT)                        :: uxy3
REAL (KIND=qPrec), INTENT(IN OUT)                        :: ux2y2

INTEGER :: n1,n2,n3,n4,m1,m2,m3,m4

INTEGER :: intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess,  &
    maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur,  &
    kcycle,iprer,ipost,intpol,kps
REAL (KIND=qPrec) :: xa,xb,yc,yd,tolmax,relmax
COMMON/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,  &
    iguess, maxcy,method,nwork,lwork,itero,ngrid,  &
    klevel,kcur,kcycle,iprer,ipost,intpol,kps
COMMON/fmud2cr/xa,xb,yc,yd,tolmax,relmax
REAL (KIND=qPrec) :: dlx,dly,dyox,dxoy,dlx2,dly2,dlxx,dlxy,dlyy,dlxy2,  &
    dlxy4,dxxxy4,dxyyy4,dxxyy,tdlx3,tdly3,dlx4,dly4, dlxxx,dlyyy
COMMON/com2dcr/dyox,dxoy,dlx2,dly2,dlxy,dlxy2,dlxy4,  &
    dxxxy4,dxyyy4,dxxyy,dlxxx,dlyyy
COMMON/pde2com/dlx,dly,dlxx,dlyy,tdlx3,tdly3,dlx4,dly4

n1=ny-1
n2=ny-2
n3=ny-3
n4=ny-4
m1=nx-1
m2=nx-2
m3=nx-3
m4=nx-4

IF (i == 1) THEN
  
  IF ((j > 2.AND.j < ny-1)) THEN
!     x=xa, yinterior
    ux3y=(5*u(1,j-1)-18*u(2,j-1)+24*u(3,j-1)-14*u(4,j-1)+3*u(5,j-1)  &
        -5*u(1,j+1)+18*u(2,j+1)-24*u(3,j+1)+14*u(4,j+1)-3*u(5,j+1)) /dxxxy4
    uxy3=(3*u(1,j-2)-4*u(2,j-2)+u(3,j-2) -6*u(1,j-1)+8*u(2,j-1)-2*u(3,j-1)  &
        +6*u(1,j+1)-8*u(2,j+1)+2*u(3,j+1) -3*u(1,j+2)+4*u(2,j+2)-u(3,j+2))/dxyyy4
  ELSE IF (j == 1) THEN
!     (xa,yc)
    ux3y=(15*u(1,1)-54*u(2,1)+72*u(3,1)-42*u(4,1)+9*u(5,1)  &
        -20*u(1,2)+72*u(2,2)-96*u(3,2)+56*u(4,2)-12*u(5,2)  &
        +5*u(1,3)-18*u(2,3)+24*u(3,3)-14*u(4,3)+3*u(5,3)) /dxxxy4
    uxy3=(15*u(1,1)-20*u(2,1)+5*u(3,1) -54*u(1,2)+72*u(2,2)-18*u(3,2)  &
        +72*u(1,3)-96*u(2,3)+24*u(3,3) -42*u(1,4)+56*u(2,4)-14*u(3,4)  &
        +9*u(1,5)-12*u(2,5)+3*u(3,5)) /dxyyy4
    ux2y2=(4*u(1,1)-10*u(2,1)+8*u(3,1)-2*u(4,1)  &
        -10*u(1,2)+25*u(2,2)-20*u(3,2)+5*u(4,2)  &
        +8*u(1,3)-20*u(2,3)+16*u(3,3)-4*u(4,3)  &
        -2*u(1,4)+5*u(2,4)-4*u(3,4)+u(4,4)) /dxxyy
  ELSE IF (j == 2) THEN
!     (xa,yc+dly)
    ux3y=(5*u(1,1)-18*u(2,1)+24*u(3,1)-14*u(4,1)+3*u(5,1)  &
        -5*u(1,3)+18*u(2,3)-24*u(3,3)+14*u(4,3)-3*u(5,3)) /dxxxy4
    uxy3=(9*u(1,1)-12*u(2,1)+3*u(3,1) -30*u(1,2)+40*u(2,2)-10*u(3,2)  &
        +36*u(1,3)-48*u(2,3)+12*u(3,3) -18*u(1,4)+24*u(2,4)-6*u(3,4)  &
        +3*u(1,5)-4*u(2,5)+u(3,5)) /dxyyy4
  ELSE IF (j == ny-1) THEN
!     x=xa,y=yd-dly
    ux3y=(5*u(1,j-1)-18*u(2,j-1)+24*u(3,j-1)-14*u(4,j-1)+3*u(5,j-1)  &
        -5*u(1,j+1)+18*u(2,j+1)-24*u(3,j+1)+14*u(4,j+1)-3*u(5,j+1))
    uxy3=(5*u(1,n2)-18*u(2,n2)+24*u(3,n2)-14*u(4,n2)+3*u(5,n2)  &
        -5*u(1,ny)+18*u(2,ny)-24*u(3,ny)+14*u(4,ny)-3*u(5,ny)) /dxyyy4
  ELSE IF (j == ny) THEN
!     x=xa, y=yd
    ux3y=(-5*u(1,n2)+18*u(2,n2)-24*u(3,n2)+14*u(4,n2)-3*u(5,n2)  &
        +20*u(1,n1)-72*u(2,n1)+96*u(3,n1)-56*u(4,n1)+12*u(5,n1)  &
        -15*u(1,ny)+54*u(2,ny)-72*u(3,ny)+42*u(4,ny)-9*u(5,ny)) /dxxxy4
    uxy3=(-9*u(1,n4)+12*u(2,n4)-3*u(3,n4) +42*u(1,n3)-56*u(2,n3)+14*u(3,n3)  &
        -72*u(1,n2)+96*u(2,n2)-24*u(3,n2) +54*u(1,n1)-72*u(2,n1)+18*u(3,n1)  &
        -15*u(1,ny)+20*u(2,ny)-5*u(3,ny)) /dxyyy4
    ux2y2=(-2*u(1,n3)+5*u(2,n3)-4*u(3,n3)+u(4,n3)  &
        +8*u(1,n2)-20*u(2,n2)+16*u(3,n2)-4*u(4,n2)  &
        -10*u(1,n1)+25*u(2,n1)-20*u(3,n1)+5*u(4,n1)  &
        +4*u(1,ny)-10*u(2,ny)+8*u(3,ny)-2*u(4,ny)) /dxxyy
  END IF
  
ELSE IF (i == 2) THEN
  
  IF ((j > 2.AND.j < ny-1)) THEN
!     x=xa+dlx, y interior
    ux3y=(3*u(1,j-1)-10*u(2,j-1)+12*u(3,j-1)-6*u(4,j-1)+u(5,j-1)  &
        -3*u(1,j+1)+10*u(2,j+1)-12*u(3,j+1)+6*u(4,j+1)-u(5,j+1))/dxxxy4
    uxy3=(u(1,j-2)-u(3,j-2)-2*u(1,j-1)+2*u(3,j-1)  &
        +2*u(1,j+1)-2*u(3,j+1)-u(1,j+2)+u(3,j+2))/dxyyy4
  ELSE IF (j == 1) THEN
!     x=xa+dlx, y=yc
    ux3y=(9*u(1,1)-30*u(2,1)+36*u(3,1)-18*u(4,1)+3*u(5,1)  &
        -12*u(1,2)+40*u(2,2)-48*u(3,2)+24*u(4,2)-4*u(5,2)  &
        +3*u(1,3)-10*u(2,3)+12*u(3,3)-6*u(4,3)+u(5,3)) /dxxxy4
    uxy3=(5*u(1,1)-5*u(3,1)-18*u(1,2)+18*u(3,2)  &
        +24*u(1,3)-24*u(3,3)-14*u(1,4) +14*u(3,4)+3*u(1,5)-3*u(3,5))  &
        /dxyyy4
  ELSE IF (j == 2) THEN
!     at x=xa+dlx,y=yc+dly
    ux3y=(3*u(1,1)-10*u(2,1)+12*u(3,1)-6*u(4,1)+u(5,1)  &
        -3*u(1,3)+10*u(2,3)-12*u(3,3)+6*u(4,3)-u(5,3)) /dxxxy4
    uxy3=(3*u(1,1)-3*u(3,1)-10*u(1,2)+10*u(3,2)  &
        +12*u(1,3)-12*u(3,3)-6*u(1,4)+6*u(3,4) +u(1,5)-u(3,5))  &
        /dxyyy4
  ELSE IF (j == ny-1) THEN
!     x=xa+dlx,y=yd-dly
    ux3y=(3*u(1,n2)-10*u(2,n2)+12*u(3,n2)-6*u(4,n2)+u(5,n2)  &
        -3*u(1,ny)+10*u(2,ny)-12*u(3,ny)+6*u(4,ny)-u(5,ny)) /dxxxy4
    uxy3=(-u(1,n4)+u(3,n4)+6*u(1,n3)-6*u(3,n3)  &
        -12*u(1,n2)+12*u(3,n2)+10*u(1,n1)-10*u(3,n1) -3*u(1,ny)+3*u(3,ny))  &
        /dxyyy4
  ELSE IF (j == ny) THEN
!     at x=xa+dlx,y=yd
    ux3y=(-3*u(1,n2)+10*u(2,n2)-12*u(3,n2)+6*u(4,n2)-u(5,n2)  &
        +12*u(1,n1)-40*u(2,n1)+48*u(3,n1)-24*u(4,n1)+4*u(5,n1)  &
        -9*u(1,ny)+30*u(2,ny)-36*u(3,ny)+18*u(4,ny)-3*u(5,ny)) /dxxxy4
    uxy3=(-3*u(1,n4)+3*u(3,n4)+14*u(1,n3)-14*u(3,n3)  &
        -24*u(1,n2)+24*u(3,n2)+18*u(1,n1)-18*u(3,n1) -5*u(1,ny)+5*u(3,ny))  &
        /dxyyy4
  END IF
  
ELSE IF (i > 2 .AND. i < nx-1) THEN
  
  IF (j == 1) THEN
!     y=yc,x interior
    ux3y=(3.0*u(i-2,1)-6.0*u(i-1,1)+6.0*u(i+1,1)-3.0*u(i+2,1)  &
        -4.0*u(i-2,2)+8.0*u(i-1,2)-8.0*u(i+1,2)+4.0*u(i+2,2)  &
        +u(i-2,3)-2.0*u(i-1,3)+2.0*u(i+1,3)-u(i+2,3)) /dxxxy4
    uxy3=(5.0*u(i-1,1)-5.0*u(i+1,1)-18.0*u(i-1,2)+18.0*u(i+1,2)  &
        +24.0*u(i-1,3)-24.0*u(i+1,3)-14.0*u(i-1,4)+14.0*u(i+1,4)  &
        +3.0*u(i-1,5)-3.0*u(i+1,5)) /dxyyy4
  ELSE IF (j == 2) THEN
!     y=yc+dly,x interior
    ux3y=(u(i-2,1)-2.0*u(i-1,1)+2.0*u(i+1,1)-u(i+2,1)  &
        -u(i-2,3)+2.0*u(i-1,3)-2.0*u(i+1,3)+u(i+2,3)) /dxxxy4
    uxy3=(u(i-1,1)-u(i+1,1)-2.0*u(i-1,2)+2.0*u(i+1,2)  &
        +2.0*u(i-1,4)-2.0*u(i+1,4)-u(i-1,5)+u(i+1,5)) /dxyyy4
  ELSE IF (j == ny-1) THEN
!     y=yd-dly, x interior
    ux3y=(u(i-2,n2)-2.0*u(i-1,n2)+2.0*u(i+1,n2)-u(i+2,n2)  &
        -u(i-2,ny)+2.0*u(i-1,ny)-2.0*u(i+1,ny)+u(i+2,ny)) /dxxxy4
    uxy3=(-u(i-1,n4)+u(i+1,n4)+6.0*u(i-1,n3)-6.0*u(i+1,n3)  &
        -12.0*u(i-1,n2)+12.0*u(i+1,n2)+10.0*u(i-1,n1)-10.0*u(i+1,n1)  &
        -3.0*u(i-1,ny)+3.0*u(i+1,ny)) /dxyyy4
  ELSE IF (j == ny) THEN
!     at y=yd, x interior
    ux3y=(-u(i-2,n2)+2.0*u(i-1,n2)-2.0*u(i+1,n2)+u(i+2,n2)  &
        +4.0*u(i-2,n1)-8.0*u(i-1,n1)+8.0*u(i+1,n1)-4.0*u(i+2,n1)  &
        -3.0*u(i-2,ny)+6.0*u(i-1,ny)-6.0*u(i+1,ny)+3.0*u(i+2,ny)) /dxxxy4
    uxy3=(-3.0*u(i-1,n4)+3.0*u(i+1,n4)+14.0*u(i-1,n3)-14.0*u(i+1,n3)  &
        -24.0*u(i-1,n2) +24.0*u(i+1,n2)+18.0*u(i-1,n1)-18.0*u(i+1,n1)  &
        -5.0*u(i-1,ny)+5.0*u(i+1,ny)) /dxyyy4
  END IF
  
ELSE IF (i == nx-1) THEN
  
  IF ((j > 2.AND.j < ny-1)) THEN
!     x=xb-dlx,y interior
    ux3y=(-u(m4,j-1)+6.*u(m3,j-1)-12.*u(m2,j-1)+10.*u(m1,j-1)-3.*u(nx  &
        ,j-1)+u(m4,j+1)-6.*u(m3,j+1)+12.*u(m2,j+1)-10.*u(m1,j+1)+3.*u(nx,j  &
        +1)) /dxxxy4
    uxy3=(u(m2,j-2)-u(nx,j-2)-2.*u(m2,j-1)+2.*u(nx,j-1)  &
        +2.*u(m2,j+1)-2.*u(nx,j+1)-u(m2,j+2)+u(nx,j+2)) /dxyyy4
  ELSE IF (j == 1) THEN
!     at x=xb-dlx, y=yc
    ux3y=(-3.0*u(m4,1)+18.0*u(m3,1)-36.0*u(m2,1)+30.0*u(m1,1)-9.0*u(  &
        nx,1)+4.0*u(m4,2)-24.0*u(m3,2)+48.0*u(m2,2)-40.0*u(m1,2)+12.0*u(nx  &
        ,2)-u(m4,3)+6.0*u(m3,3)-12.0*u(m2,3)+10.0*u(m1,3)-3.0*u(nx,3)) /dxxxy4
    uxy3=(5.0*u(m2,1)-5.0*u(nx,1)-18.0*u(m2,2)+18.0*u(nx,2)  &
        +24.0*u(m2,3)-24.0*u(nx,3)-14.0*u(m2,4)+14.0*u(nx,4)  &
        +3.0*u(m2,5)-3.0*u(nx,5)) /dxyyy4
  ELSE IF (j == 2) THEN
!     x=xb-dlx,y=yc+dly
    ux3y=(-u(m4,1)+6.0*u(m3,1)-12.0*u(m2,1)+10.*u(m1,1)-3.*u(nx,1)  &
        +u(m4,3)-6.0*u(m3,3)+12.0*u(m2,3)-10.*u(m1,3)+3.*u(nx,3)) /dxxxy4
    uxy3=(3.0*u(m2,1)-3.*u(nx,1)-10.*u(m2,2)+10.*u(nx,2)  &
        +12.*u(m2,3)-12.*u(nx,3)-6.*u(m2,4)+6.*u(nx,4) +u(m2,5)-u(nx,5)) / dxyyy4
  ELSE IF (j == ny-1) THEN
!     at x=xb-dlx,y=yd-dly
    ux3y=(-u(m4,n2)+6.*u(m3,n2)-12.*u(m2,n2)+10.*u(m1,n2)-3.*u(nx,n2)  &
        +u(m4,ny)-6.*u(m3,ny)+12.*u(m2,ny)-10.*u(m1,ny)+3.*u(nx,ny)) /dxxxy4
    uxy3=(-u(m2,n4)+u(nx,n4)+6*u(m2,n3)-6.*u(nx,n3)  &
        -12.*u(m2,n2)+12.*u(nx,n2)+10.*u(m2,n1)-10.*u(nx,n1)  &
        -3.*u(m2,ny)+3.*u(nx,ny)) / dxyyy4
  ELSE IF (j == ny) THEN
!     at x=xb.dlx,y=yd
    ux3y=(u(m4,n2)-6.*u(m3,n2)+12.*u(m2,n2)-10.*u(m1,n2)+3.*u(nx,n2)  &
        -4.*u(m4,n1)+24.*u(m3,n1)-48.*u(m2,n1)+40.*u(m1,n1)-12.*u(nx,n1)  &
        +3.*u(m4,ny)-18.*u(m3,ny)+36.*u(m2,ny)-30.*u(m1,ny)+9.*u(nx,ny)) / dxxxy4
    uxy3=(-3.*u(m2,n4)+3.*u(nx,n4)+14.*u(m2,n3)-14.*u(nx,n3)  &
        -24.*u(m2,n2)+24.*u(nx,n2)+18.*u(m2,n1)-18.*u(nx,n1)  &
        -5.*u(m2,ny)+5.*u(nx,ny)) / dxyyy4
  END IF
  
ELSE IF (i == nx) THEN
  
  IF ((j > 2.AND.j < ny-1)) THEN
!     x=xb,y interior
    ux3y=(-3.*u(m4,j-1)+14.*u(m3,j-1)-24.*u(m2,j-1)+18.*u(m1,j-1)-5.*  &
        u(nx,j-1)+3.*u(m4,j+1)-14.*u(m3,j+1)+24.*u(m2,j+1)-18.*u(m1,j+1)+  &
        5.*u(nx,j+1)) / dxxxy4
    uxy3=(-u(m2,j-2)+4.*u(m1,j-2)-3.*u(nx,j-2)  &
        +2.*u(m2,j-1)-8.*u(m1,j-1)+6.*u(nx,j-1)  &
        -2.*u(m2,j+1)+8.*u(m1,j+1)-6.*u(nx,j+1)  &
        +u(m2,j+2)-4.*u(m1,j+2)+3.*u(nx,j+2)) / dxyyy4
  ELSE IF (j == 1) THEN
!     x=xb,y=yc
    ux3y=(-9.*u(m4,1)+42.*u(m3,1)-72.*u(m2,1)+54.*u(m1,1)-15.*u(nx,1)  &
        +12.*u(m4,2)-56.*u(m3,2)+96.*u(m2,2)-72.*u(m1,2)+20.*u(nx,2)  &
        -3.*u(m4,3)+14.*u(m3,3)-24.*u(m2,3)+18.*u(m1,3)-5.*u(nx,3)) /dxxxy4
    uxy3=(-5.*u(m2,1)+20.*u(m1,1)-15.*u(nx,1)  &
        +18.*u(m2,2)-72.*u(m1,2)+54.*u(nx,2)  &
        -24.*u(m2,3)+96.*u(m1,3)-72.*u(nx,3)  &
        +14.*u(m2,4)-56.*u(m1,4)+42.*u(nx,4)  &
        -3.*u(m2,5)+12.*u(m1,5)-9.*u(nx,5)) / dxyyy4
    ux2y2=(-2.*u(m3,1)+8.*u(m2,1)-10.*u(m1,1)+4.*u(nx,1)  &
        +5.*u(m3,2)-20.*u(m2,2)+25.*u(m1,2)-10.*u(nx,2)  &
        -4.*u(m3,3)+16.*u(m2,3)-20.*u(m1,3)+8.*u(nx,3)  &
        +u(m3,4)-4.*u(m2,4)+5.*u(m1,4)-2.*u(nx,4)) / dxxyy
  ELSE IF (j == 2) THEN
!     x=xb,y=yc+dly
    ux3y=(-3.*u(m4,1)+14.*u(m3,1)-24.*u(m2,1)+18.*u(m1,1)-5.*u(nx,1)  &
        +3.*u(m4,3)-14.*u(m3,3)+24.*u(m2,3)-18.*u(m1,3)+5.*u(nx,3)) / dxxxy4
    uxy3=(-3.*u(m2,1)+12.*u(m1,1)-9.*u(nx,1)  &
        +10.*u(m2,2)-40.*u(m1,2)+30.*u(nx,2)  &
        -12.*u(m2,3)+48.*u(m1,3)-36.*u(nx,3)  &
        +6.*u(m2,4)-24.*u(m1,4)+18.*u(nx,4)  &
        -u(m2,5)+4.*u(m1,5)-3.*u(nx,5)) / dxyyy4
  ELSE IF (j == ny-1) THEN
!     x=xb,y=yd-dly
    ux3y=(-3.*u(m4,n2)+14.*u(m3,n2)-24.*u(m2,n2)+18.*u(m1,n2)-5.*u(nx  &
        ,n2)+3.*u(m4,ny)-14.*u(m3,ny)+24.*u(m2,ny)-18.*u(m1,ny)+5.*u(nx,ny  &
        )) / dxxxy4
    uxy3=(u(m2,n4)-4.*u(m1,n4)+3.*u(nx,n4)  &
        -6.*u(m2,n3)+24.*u(m1,n3)-18.*u(nx,n3)  &
        +12.*u(m2,n2)-48.*u(m1,n2)+36.*u(nx,n2)  &
        -10.*u(m2,n1)+40.*u(m1,n1)-30.*u(nx,n1)  &
        +3.*u(m2,ny)-12.*u(m1,ny)+9.*u(nx,ny)) / dxyyy4
  ELSE IF (j == ny) THEN
!     x=xb,y=yd
    ux3y=(3.*u(m4,n2)-14.*u(m3,n2)+24.*u(m2,n2)-18.*u(m1,n2)+5.*u(nx,  &
        n2)-12.*u(m4,n1)+56.*u(m3,n1)-96.*u(m2,n1)+72.*u(m1,n1)-20.*u(nx,  &
        n1)+9.*u(m4,ny)-42.*u(m3,ny)+72.*u(m2,ny)-54.*u(m1,ny)+15.*u(nx,ny  &
        )) / dxxxy4
    uxy3=(3.*u(m2,n4)-12.*u(m1,n4)+9.*u(nx,n4)  &
        -14.*u(m2,n3)+56.*u(m1,n3)-42.*u(nx,n3)  &
        +24.*u(m2,n2)-96.*u(m1,n2)+72.*u(nx,n2)  &
        -18.*u(m2,n1)+72.*u(m1,n1)-54.*u(nx,n1)  &
        +5.*u(m2,ny)-20.*u(m1,ny)+15.*u(nx,ny)) / dxyyy4
    ux2y2=(u(m3,n3)-4.*u(m2,n3)+5.*u(m1,n3)-2.*u(nx,n3)  &
        -4.*u(m3,n2)+16.*u(m2,n2)-20.*u(m1,n2)+8.*u(nx,n2)  &
        +5.0*u(m3,n1)-20.*u(m2,n1)+25.*u(m1,n1)-10.*u(nx,n1)  &
        -2.*u(m3,ny)+8.*u(m2,ny)-10.*u(m1,ny)+4.*u(nx,ny)) / dxxyy
  END IF
  
END IF

RETURN
END SUBROUTINE pde2cr

SUBROUTINE pde3(nx,ny,nz,u,i,j,k,ux3,ux4,uy3,uy4,uz3,uz4, nxa,nyc,nze)

!     estimate third and fourth partial derivatives in x,y,z

IMPLICIT NONE
INTEGER, INTENT(IN OUT)                  :: nx
INTEGER, INTENT(IN OUT)                  :: ny
INTEGER, INTENT(IN OUT)                      :: nz
REAL (KIND=qPrec), INTENT(IN OUT)                         :: u(nx,ny,nz)
INTEGER, INTENT(IN OUT)                  :: i
INTEGER, INTENT(IN OUT)                      :: j
INTEGER, INTENT(IN OUT)                      :: k
REAL (KIND=qPrec), INTENT(IN OUT)                     :: ux3
REAL (KIND=qPrec), INTENT(IN OUT)                     :: ux4
REAL (KIND=qPrec), INTENT(IN OUT)                     :: uy3
REAL (KIND=qPrec), INTENT(IN OUT)                     :: uy4
REAL (KIND=qPrec), INTENT(IN OUT)                        :: uz3
REAL (KIND=qPrec), INTENT(IN OUT)                        :: uz4
INTEGER, INTENT(IN OUT)                  :: nxa
INTEGER, INTENT(IN OUT)                  :: nyc
INTEGER, INTENT(IN OUT)                  :: nze



REAL (KIND=qPrec) :: dlx,dly,dlz,dlxx,dlyy,dlzz,tdlx3,tdly3,tdlz3,dlx4,dly4,dlz4
COMMON/pde3com/dlx,dly,dlz,dlxx,dlyy,dlzz,tdlx3,tdly3,tdlz3, dlx4,dly4,dlz4


!     x,y partial derivatives

CALL p3de2(nx,ny,u(1,1,k),i,j,ux3,ux4,uy3,uy4,nxa,nyc)

!     z partial derivatives

IF (nze /= 0) THEN
  
!     nonperiodic in z
  
  IF(k > 2 .AND. k < nz-1) THEN
    uz3=(-u(i,j,k-2)+2.0*u(i,j,k-1)-2.0*u(i,j,k+1)+u(i,j,k+2))/tdlz3
    uz4=(u(i,j,k-2)-4.0*u(i,j,k-1)+6.0*u(i,j,k)-4.0*u(i,j,k+1)+  &
        u(i,j,k+2))/dlz4
  ELSE IF (k == 1) THEN
    uz3=(-5.0*u(i,j,1)+18.0*u(i,j,2)-24.0*u(i,j,3)+14.0*u(i,j,4)-  &
        3.0*u(i,j,5))/tdlz3
    uz4 = (3.0*u(i,j,1)-14.0*u(i,j,2)+26.0*u(i,j,3)-24.0*u(i,j,4)+  &
        11.0*u(i,j,5)-2.0*u(i,j,6))/dlz4
  ELSE IF (k == 2) THEN
    uz3 = (-3.0*u(i,j,1)+10.0*u(i,j,2)-12.0*u(i,j,3)+6.0*u(i,j,4)-  &
        u(i,j,5))/tdlz3
    uz4 = (2.0*u(i,j,1)-9.0*u(i,j,2)+16.0*u(i,j,3)-14.0*u(i,j,4)+6.0*  &
        u(i,j,5)-u(i,j,6))/dlz4
  ELSE IF (k == nz-1) THEN
    uz3 = (u(i,j,nz-4)-6.0*u(i,j,nz-3)+12.0*u(i,j,nz-2)-10.0*  &
        u(i,j,nz-1)+3.0*u(i,j,nz))/tdlz3
    uz4 = (-u(i,j,nz-5)+6.0*u(i,j,nz-4)-14.0*u(i,j,nz-3)+16.0*  &
        u(i,j,nz-2)-9.0*u(i,j,nz-1)+2.0*u(i,j,nz))/dlz4
  ELSE IF (k == nz) THEN
    uz3 = (3.0*u(i,j,nz-4)-14.0*u(i,j,nz-3)+24.0*u(i,j,nz-2)-18.0*  &
        u(i,j,nz-1)+5.0*u(i,j,nz))/tdlz3
    uz4 = (-2.0*u(i,j,nz-5)+11.0*u(i,j,nz-4)-24.0*u(i,j,nz-3)+26.0*  &
        u(i,j,nz-2)-14.0*u(i,j,nz-1)+3.0*u(i,j,nz))/dlz4
  END IF
ELSE
  
!     periodic in z so use symmetric formula even "near" z boundaies
  
  IF(k > 2 .AND. k < nz-1) THEN
    uz3=(-u(i,j,k-2)+2.0*u(i,j,k-1)-2.0*u(i,j,k+1)+u(i,j,k+2))/tdlz3
    uz4=(u(i,j,k-2)-4.0*u(i,j,k-1)+6.0*u(i,j,k)-4.0*u(i,j,k+1)+  &
        u(i,j,k+2))/dlz4
  ELSE IF (k == 1) THEN
    uz3 = (-u(i,j,nz-2)+2.0*u(i,j,nz-1)-2.0*u(i,j,2)+u(i,j,3))/tdlz3
    uz4 = (u(i,j,nz-2)-4.0*u(i,j,nz-1)+6.0*u(i,j,1)-4.0*u(i,j,2)+  &
        u(i,j,3))/dlz4
  ELSE IF (k == 2) THEN
    uz3 = (-u(i,j,nz-1)+2.0*u(i,j,1)-2.0*u(i,j,3)+u(i,j,4))/(tdlz3)
    uz4 = (u(i,j,nz-1)-4.0*u(i,j,1)+6.0*u(i,j,2)-4.0*u(i,j,3)+ u(i,j,4))/dlz4
  ELSE IF (k == nz-1) THEN
    uz3 = (-u(i,j,nz-3)+2.0*u(i,j,nz-2)-2.0*u(i,j,1)+u(i,j,2))/tdlz3
    uz4 = (u(i,j,nz-3)-4.0*u(i,j,nz-2)+6.0*u(i,j,nz-1)-4.0*u(i,j,1)+  &
        u(i,j,2))/ dlz4
  ELSE IF (k == nz) THEN
    uz3 = (-u(i,j,nz-2)+2.0*u(i,j,nz-1)-2.0*u(i,j,2)+u(i,j,3))/tdlz3
    uz4 = (u(i,j,nz-2)-4.0*u(i,j,nz-1)+6.0*u(i,j,nz)-4.0*u(i,j,2)+  &
        u(i,j,3))/dlz4
  END IF
END IF
RETURN
END SUBROUTINE pde3

SUBROUTINE p3de2(nx,ny,u,i,j,ux3,ux4,uy3,uy4,nxa,nyc)

!     third and fourth partial derivatives in x and y

IMPLICIT NONE
INTEGER, INTENT(IN OUT)                  :: nx
INTEGER, INTENT(IN OUT)                      :: ny
REAL (KIND=qPrec), INTENT(IN OUT)                         :: u(nx,ny)
INTEGER, INTENT(IN OUT)                  :: i
INTEGER, INTENT(IN OUT)                      :: j
REAL (KIND=qPrec), INTENT(IN OUT)                     :: ux3
REAL (KIND=qPrec), INTENT(IN OUT)                     :: ux4
REAL (KIND=qPrec), INTENT(IN OUT)                        :: uy3
REAL (KIND=qPrec), INTENT(IN OUT)                        :: uy4
INTEGER, INTENT(IN OUT)                  :: nxa
INTEGER, INTENT(IN OUT)                  :: nyc

INTEGER :: l

REAL (KIND=qPrec) :: dlx,dly,dlz,dlxx,dlyy,dlzz,tdlx3,tdly3,tdlz3,dlx4,dly4,dlz4
COMMON/pde3com/dlx,dly,dlz,dlxx,dlyy,dlzz,tdlx3,tdly3,tdlz3, dlx4,dly4,dlz4


l=ny

!     x partial derivatives

CALL p3de1(nx,u(1,j),i,ux3,ux4,nxa)

!     y partial derivatives

IF (nyc /= 0) THEN
  
!     not periodic in y
  
  IF (j > 2 .AND. j < ny-1) THEN
    uy3 = (-u(i,j-2)+2.0*u(i,j-1)-2.0*u(i,j+1)+u(i,j+2))/tdly3
    uy4 = (u(i,j-2)-4.0*u(i,j-1)+6.0*u(i,j)-4.0*u(i,j+1)+u(i,j+2))/ dly4
  ELSE IF (j == 1) THEN
    uy3 = (-5.0*u(i,1)+18.0*u(i,2)-24.0*u(i,3)+14.0*u(i,4)- 3.0*u(i,5))/tdly3
    uy4 = (3.0*u(i,1)-14.0*u(i,2)+26.0*u(i,3)-24.0*u(i,4)+  &
        11.0*u(i,5)-2.0*u(i,6))/dly4
  ELSE IF (j == 2) THEN
    uy3 = (-3.0*u(i,1)+10.0*u(i,2)-12.0*u(i,3)+6.0*u(i,4)-u(i,5))/ tdly3
    uy4 = (2.0*u(i,1)-9.0*u(i,2)+16.0*u(i,3)-14.0*u(i,4)+6.0*u(i,5)-  &
        u(i,6))/dly4
  ELSE IF (j == ny-1) THEN
    uy3 = (u(i,l-4)-6.0*u(i,l-3)+12.0*u(i,l-2)-10.0*u(i,l-1)+  &
        3.0*u(i,l))/tdly3
    uy4 = (-u(i,l-5)+6.0*u(i,l-4)-14.0*u(i,l-3)+16.0*u(i,l-2)-  &
        9.0*u(i,l-1)+2.0*u(i,l))/dly4
  ELSE IF (j == ny) THEN
    uy3 = (3.0*u(i,l-4)-14.0*u(i,l-3)+24.0*u(i,l-2)-18.0*u(i,l-1)+  &
        5.0*u(i,l))/tdly3
    uy4 = (-2.0*u(i,l-5)+11.0*u(i,l-4)-24.0*u(i,l-3)+26.0*u(i,l-2)-  &
        14.0*u(i,l-1)+3.0*u(i,l))/dly4
  END IF
ELSE
  
!     periodic in y
  
  IF (j > 2 .AND. j < ny-1) THEN
    uy3 = (-u(i,j-2)+2.0*u(i,j-1)-2.0*u(i,j+1)+u(i,j+2))/tdly3
    uy4 = (u(i,j-2)-4.0*u(i,j-1)+6.0*u(i,j)-4.0*u(i,j+1)+u(i,j+2))/ dly4
  ELSE IF (j == 1) THEN
    uy3 = (-u(i,l-2)+2.0*u(i,l-1)-2.0*u(i,2)+u(i,3))/tdly3
    uy4 = (u(i,l-2)-4.0*u(i,l-1)+6.0*u(i,1)-4.0*u(i,2)+u(i,3))/dly4
  ELSE IF (j == 2) THEN
    uy3 = (-u(i,l-1)+2.0*u(i,1)-2.0*u(i,3)+u(i,4))/(tdly3)
    uy4 = (u(i,l-1)-4.0*u(i,1)+6.0*u(i,2)-4.0*u(i,3)+u(i,4))/dly4
  ELSE IF (j == ny-1) THEN
    uy3 = (-u(i,l-3)+2.0*u(i,l-2)-2.0*u(i,1)+u(i,2))/tdly3
    uy4 = (u(i,l-3)-4.0*u(i,l-2)+6.0*u(i,l-1)-4.0*u(i,1)+u(i,2))/ dly4
  ELSE IF (j == ny) THEN
    uy3 = (-u(i,l-2)+2.0*u(i,l-1)-2.0*u(i,2)+u(i,3))/tdly3
    uy4 = (u(i,l-2)-4.0*u(i,l-1)+6.0*u(i,l)-4.0*u(i,2)+u(i,3))/dly4
  END IF
END IF
RETURN
END SUBROUTINE p3de2

SUBROUTINE p3de1(nx,u,i,ux3,ux4,nxa)

!     third and fourth derivatives in x

IMPLICIT NONE
INTEGER, INTENT(IN OUT)                      :: nx
REAL (KIND=qPrec), INTENT(IN OUT)                         :: u(nx)
INTEGER, INTENT(IN OUT)                      :: i
REAL (KIND=qPrec), INTENT(IN OUT)                        :: ux3
REAL (KIND=qPrec), INTENT(IN OUT)                        :: ux4
INTEGER, INTENT(IN OUT)                  :: nxa

INTEGER :: k

REAL (KIND=qPrec) :: dlx,dly,dlz,dlxx,dlyy,dlzz,tdlx3,tdly3,tdlz3,dlx4,dly4,dlz4
COMMON/pde3com/dlx,dly,dlz,dlxx,dlyy,dlzz,tdlx3,tdly3,tdlz3, dlx4,dly4,dlz4


k = nx
IF (nxa /= 0) THEN
  
!     nonperiodic in x
  
  IF(i > 2 .AND. i < nx-1) THEN
    ux3 = (-u(i-2)+2.0*u(i-1)-2.0*u(i+1)+u(i+2))/tdlx3
    ux4 = (u(i-2)-4.0*u(i-1)+6.0*u(i)-4.0*u(i+1)+u(i+2))/dlx4
  ELSE IF (i == 1) THEN
    ux3 = (-5.0*u(1)+18.0*u(2)-24.0*u(3)+14.0*u(4)-3.0*u(5))/tdlx3
    ux4 = (3.0*u(1)-14.0*u(2)+26.0*u(3)-24.0*u(4)+11.0*u(5)-2.0*u(6)) /dlx4
  ELSE IF (i == 2) THEN
    ux3 = (-3.0*u(1)+10.0*u(2)-12.0*u(3)+6.0*u(4)-u(5))/tdlx3
    ux4 = (2.0*u(1)-9.0*u(2)+16.0*u(3)-14.0*u(4)+6.0*u(5)-u(6))/dlx4
  ELSE IF (i == nx-1) THEN
    ux3 = (u(k-4)-6.0*u(k-3)+12.0*u(k-2)-10.0*u(k-1)+3.0*u(k))/tdlx3
    ux4 = (-u(k-5)+6.0*u(k-4)-14.0*u(k-3)+16.0*u(k-2)-9.0*u(k-1)+  &
        2.0*u(k))/dlx4
  ELSE IF (i == nx) THEN
    ux3 = (3.0*u(k-4)-14.0*u(k-3)+24.0*u(k-2)-18.0*u(k-1)+5.0*u(k))/ tdlx3
    ux4 = (-2.0*u(k-5)+11.0*u(k-4)-24.0*u(k-3)+26.0*u(k-2)-  &
        14.0*u(k-1)+3.0*u(k))/dlx4
  END IF
ELSE
  
!     periodic in x
  
  IF(i > 2 .AND. i < nx-1) THEN
    ux3 = (-u(i-2)+2.0*u(i-1)-2.0*u(i+1)+u(i+2))/tdlx3
    ux4 = (u(i-2)-4.0*u(i-1)+6.0*u(i)-4.0*u(i+1)+u(i+2))/dlx4
  ELSE IF (i == 1) THEN
    ux3 = (-u(k-2)+2.0*u(k-1)-2.0*u(2)+u(3))/tdlx3
    ux4 = (u(k-2)-4.0*u(k-1)+6.0*u(1)-4.0*u(2)+u(3))/dlx4
  ELSE IF (i == 2) THEN
    ux3 = (-u(k-1)+2.0*u(1)-2.0*u(3)+u(4))/(tdlx3)
    ux4 = (u(k-1)-4.0*u(1)+6.0*u(2)-4.0*u(3)+u(4))/dlx4
  ELSE IF (i == nx-1) THEN
    ux3 = (-u(k-3)+2.0*u(k-2)-2.0*u(1)+u(2))/tdlx3
    ux4 = (u(k-3)-4.0*u(k-2)+6.0*u(k-1)-4.0*u(1)+u(2))/dlx4
  ELSE IF (i == nx) THEN
    ux3 = (-u(k-2)+2.0*u(k-1)-2.0*u(2)+u(3))/tdlx3
    ux4 = (u(k-2)-4.0*u(k-1)+6.0*u(k)-4.0*u(2)+u(3))/dlx4
  END IF
END IF
RETURN
END SUBROUTINE p3de1


!     factri and factrip are:
!     subroutines called by any real mudpack solver which uses line
!     relaxation(s) within multigrid iteration.  these subroutines do
!     a vectorized factorization of m simultaneous tridiagonal systems
!     of order n arising from nonperiodic or periodic discretizations

SUBROUTINE factri(m,n,a,b,c)

!     factor the m simultaneous tridiagonal systems of order n

IMPLICIT NONE
INTEGER, INTENT(IN)                      :: m
INTEGER, INTENT(IN)                      :: n
REAL (KIND=qPrec), INTENT(IN OUT)                        :: a(n,m)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: b(n,m)
REAL (KIND=qPrec), INTENT(IN OUT)                         :: c(n,m)

INTEGER :: i,j


DO i=2,n
  DO j=1,m
    a(i-1,j) = a(i-1,j)/b(i-1,j)
    b(i,j) = b(i,j)-a(i-1,j)*c(i-1,j)
  END DO
END DO
RETURN
END SUBROUTINE factri



SUBROUTINE factrp(m,n,a,b,c,d,e,sum)

!     factor the m simultaneous "tridiagonal" systems of order n
!     from discretized periodic system (leave out periodic n point)
!     (so sweeps below only go from i=1,2,...,n-1) n > 3 is necessary

IMPLICIT NONE
INTEGER, INTENT(IN)                      :: m
INTEGER, INTENT(IN)                      :: n
REAL (KIND=qPrec), INTENT(IN OUT)                     :: a(n,m)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: b(n,m)
REAL (KIND=qPrec), INTENT(IN OUT)                         :: c(n,m)
REAL (KIND=qPrec), INTENT(IN OUT)                        :: d(n,m)
REAL (KIND=qPrec), INTENT(IN OUT)                        :: e(n,m)
REAL (KIND=qPrec), INTENT(IN OUT)                        :: sum(m)

INTEGER :: i,j


DO j=1,m
  d(1,j) = a(1,j)
END DO
DO i=2,n-2
  DO j=1,m
    a(i,j) = a(i,j)/b(i-1,j)
    b(i,j) = b(i,j)-a(i,j)*c(i-1,j)
    d(i,j) = -a(i,j)*d(i-1,j)
  END DO
END DO

!     correct computation of last d element

DO j=1,m
  d(n-2,j) = c(n-2,j)+d(n-2,j)
END DO
DO j=1,m
  e(1,j) = c(n-1,j)/b(1,j)
END DO
DO i=2,n-3
  DO j=1,m
    e(i,j) = -e(i-1,j)*c(i-1,j)/b(i,j)
  END DO
END DO
DO j=1,m
  e(n-2,j) = (a(n-1,j)-e(n-3,j)*c(n-3,j))/b(n-2,j)
END DO

!     compute  inner product (e,d) for each j in sum(j)

DO j=1,m
  sum(j) = 0.
END DO
DO i=1,n-2
  DO j=1,m
    sum(j) = sum(j)+e(i,j)*d(i,j)
  END DO
END DO

!     set last diagonal element

DO j=1,m
  b(n-1,j) = b(n-1,j)-sum(j)
END DO
RETURN
END SUBROUTINE factrp

SUBROUTINE transp(n,amat)

!     transpose n by n real matrix

IMPLICIT NONE
INTEGER, INTENT(IN OUT)                      :: n
REAL (KIND=qPrec), INTENT(IN OUT)                     :: amat(n,n)

INTEGER :: i,j
REAL (KIND=qPrec) :: temp

DO i=1,n-1
  DO j=i+1,n
    temp = amat(i,j)
    amat(i,j) = amat(j,i)
    amat(j,i) = temp
  END DO
END DO
RETURN
END SUBROUTINE transp

SUBROUTINE sgfa (a,lda,n,ipvt,info)

REAL (KIND=qPrec), INTENT(IN OUT)                     :: a(lda,1)
INTEGER, INTENT(IN OUT)                  :: lda
INTEGER, INTENT(IN OUT)                      :: n
INTEGER, INTENT(IN OUT)                     :: ipvt(1)
INTEGER, INTENT(IN OUT)                     :: info


REAL (KIND=qPrec) :: t
INTEGER :: isfmax,j,k,kp1,l,nm1

info = 0
nm1 = n - 1
IF (nm1 < 1) GO TO 70
DO  k = 1, nm1
  kp1 = k + 1
  l = isfmax(n-k+1,a(k,k),1) + k - 1
  ipvt(k) = l
  IF (a(l,k) == 0.0E0) GO TO 40
  IF (l == k) GO TO 10
  t = a(l,k)
  a(l,k) = a(k,k)
  a(k,k) = t
  10       CONTINUE
  t = -1.0E0/a(k,k)
  CALL sscl(n-k,t,a(k+1,k),1)
  DO  j = kp1, n
    t = a(l,j)
    IF (l == k) GO TO 20
    a(l,j) = a(k,j)
    a(k,j) = t
    20          CONTINUE
    CALL sxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
  END DO
  GO TO 50
  40    CONTINUE
  info = k
  50    CONTINUE
END DO
70 CONTINUE
ipvt(n) = n
IF (a(n,n) == 0.0E0) info = n
RETURN
END SUBROUTINE sgfa

SUBROUTINE sgsl (a,lda,n,ipvt,b,job)

REAL (KIND=qPrec), INTENT(IN OUT)                         :: a(lda,1)
INTEGER, INTENT(IN OUT)                  :: lda
INTEGER, INTENT(IN OUT)                      :: n
INTEGER, INTENT(IN OUT)                      :: ipvt(1)
REAL (KIND=qPrec), INTENT(IN OUT)                     :: b(1)
INTEGER, INTENT(IN OUT)                  :: job


REAL (KIND=qPrec) :: sdt,t
INTEGER :: k,kb,l,nm1

nm1 = n - 1
IF (job /= 0) GO TO 50
IF (nm1 < 1) GO TO 30
DO  k = 1, nm1
  l = ipvt(k)
  t = b(l)
  IF (l == k) GO TO 10
  b(l) = b(k)
  b(k) = t
  10       CONTINUE
  CALL sxpy(n-k,t,a(k+1,k),1,b(k+1),1)
END DO
30    CONTINUE
DO  kb = 1, n
  k = n + 1 - kb
  b(k) = b(k)/a(k,k)
  t = -b(k)
  CALL sxpy(k-1,t,a(1,k),1,b(1),1)
END DO
GO TO 100
50 CONTINUE
DO  k = 1, n
  t = sdt(k-1,a(1,k),1,b(1),1)
  b(k) = (b(k) - t)/a(k,k)
END DO
IF (nm1 < 1) GO TO 90
DO  kb = 1, nm1
  k = n - kb
  b(k) = b(k) + sdt(n-k,a(k+1,k),1,b(k+1),1)
  l = ipvt(k)
  IF (l == k) GO TO 70
  t = b(l)
  b(l) = b(k)
  b(k) = t
  70       CONTINUE
END DO
90    CONTINUE
100 CONTINUE
RETURN
END SUBROUTINE sgsl

REAL (KIND=qPrec) FUNCTION sdt(n,sx,incx,sy,incy)

INTEGER, INTENT(IN )                      :: n
REAL (KIND=qPrec), INTENT(IN OUT)                         :: sx(1)
INTEGER, INTENT(IN )                      :: incx
REAL (KIND=qPrec), INTENT(IN OUT)                         :: sy(1)
INTEGER, INTENT(IN )                      :: incy
REAL (KIND=qPrec) :: stemp
INTEGER :: i, ix,iy,m,mp1

stemp = 0.0E0
sdt = 0.0E0
IF(n <= 0)RETURN
IF(incx == 1.AND.incy == 1)GO TO 20
ix = 1
iy = 1
IF(incx < 0)ix = (-n+1)*incx + 1
IF(incy < 0)iy = (-n+1)*incy + 1
DO  i = 1,n
  stemp = stemp + sx(ix)*sy(iy)
  ix = ix + incx
  iy = iy + incy
END DO
sdt = stemp
RETURN
20 m = MOD(n,5)
IF( m == 0 ) GO TO 40
DO  i = 1,m
  stemp = stemp + sx(i)*sy(i)
END DO
IF( n < 5 ) GO TO 60
40 mp1 = m + 1
DO  i = mp1,n,5
  stemp = stemp + sx(i)*sy(i) + sx(i + 1)*sy(i + 1) +  &
      sx(i + 2)*sy(i + 2) + sx(i + 3)*sy(i + 3) + sx(i + 4)*sy(i + 4)
END DO
60 sdt = stemp
RETURN
END FUNCTION sdt

INTEGER FUNCTION isfmax(n,sx,incx)

INTEGER, INTENT(IN )                      :: n
REAL (KIND=qPrec), INTENT(IN OUT)                     :: sx(1)
INTEGER, INTENT(IN )                      :: incx
REAL (KIND=qPrec) :: smax
INTEGER :: i, ix

isfmax = 0
IF( n < 1 ) RETURN
isfmax = 1
IF(n == 1)RETURN
IF(incx == 1)GO TO 20
ix = 1
smax = ABS(sx(1))
ix = ix + incx
DO  i = 2,n
  IF(ABS(sx(ix)) <= smax) GO TO 5
  isfmax = i
  smax = ABS(sx(ix))
  5    ix = ix + incx
END DO
RETURN
20 smax = ABS(sx(1))
DO  i = 2,n
  IF(ABS(sx(i)) <= smax) CYCLE
  isfmax = i
  smax = ABS(sx(i))
END DO
RETURN
END FUNCTION isfmax

SUBROUTINE sxpy(n,sa,sx,incx,sy,incy)

INTEGER, INTENT(IN )                      :: n
REAL (KIND=qPrec), INTENT(IN OUT)                         :: sa
REAL (KIND=qPrec), INTENT(IN OUT)                         :: sx(1)
INTEGER, INTENT(IN )                      :: incx
REAL (KIND=qPrec), INTENT(IN OUT)                        :: sy(1)
INTEGER, INTENT(IN )                      :: incy

INTEGER :: i, ix,iy,m,mp1

IF(n <= 0)RETURN
IF (sa == 0.0) RETURN
IF(incx == 1.AND.incy == 1)GO TO 20
ix = 1
iy = 1
IF(incx < 0)ix = (-n+1)*incx + 1
IF(incy < 0)iy = (-n+1)*incy + 1
DO  i = 1,n
  sy(iy) = sy(iy) + sa*sx(ix)
  ix = ix + incx
  iy = iy + incy
END DO
RETURN
20 m = MOD(n,4)
IF( m == 0 ) GO TO 40
DO  i = 1,m
  sy(i) = sy(i) + sa*sx(i)
END DO
IF( n < 4 ) RETURN
40 mp1 = m + 1
DO  i = mp1,n,4
  sy(i) = sy(i) + sa*sx(i)
  sy(i + 1) = sy(i + 1) + sa*sx(i + 1)
  sy(i + 2) = sy(i + 2) + sa*sx(i + 2)
  sy(i + 3) = sy(i + 3) + sa*sx(i + 3)
END DO
RETURN
END SUBROUTINE sxpy

SUBROUTINE sscl(n,sa,sx,incx)

INTEGER, INTENT(IN )                      :: n
REAL (KIND=qPrec), INTENT(IN OUT)                         :: sa
REAL (KIND=qPrec), INTENT(IN OUT)                        :: sx(1)
INTEGER, INTENT(IN )                      :: incx

INTEGER :: i, m,mp1, nincx

IF(n <= 0)RETURN
IF(incx == 1)GO TO 20
nincx = n*incx
DO  i = 1,nincx,incx
  sx(i) = sa*sx(i)
END DO
RETURN
20 m = MOD(n,5)
IF( m == 0 ) GO TO 40
DO  i = 1,m
  sx(i) = sa*sx(i)
END DO
IF( n < 5 ) RETURN
40 mp1 = m + 1
DO  i = mp1,n,5
  sx(i) = sa*sx(i)
  sx(i + 1) = sa*sx(i + 1)
  sx(i + 2) = sa*sx(i + 2)
  sx(i + 3) = sa*sx(i + 3)
  sx(i + 4) = sa*sx(i + 4)
END DO
RETURN
END SUBROUTINE sscl



END MODULE libinterp
